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1.  INTROtflrtTFION  AND  SUMMARY 

■  -  -  i'  i  ^  X  Oil 

This  final  report  summarizes  work  performed  on  the  Adaptive  Decentralized 
Control  project  (under  contract  F4920-81-C-0051)  during  the  period  June  1981  - 
July  1984.  The  objective  of  this  research  effort  was  the  development  of  a  new 
concept  for  the  design  of  decentralized  controllers  for  large  scale  systems. 

The  modeling,  analysis  and  control  of  large-scale  systems  Is  an 
Increasingly  Important  problem  In  such  diverse  areas  as  defense  systems, 
communication  and  computer  networks  and  transportatl on  systems.  The  size  and 
complexity  of  many  systems  make  It  difficult  or  Impractical  to  use  centralized 
control  structures.  Furthermore,  considerations  of  communication  costs, 
system  reliability,  computational  requirements  and  response  time  provide 
strong  Incentives  for  the  use  of  distributed  control  architectures.  The  basic 
focus  of  our  research  Is  on  a  framework  within  which  decentralized  controller 
structures  can  be  analyzed  and  developed.  The  motivation  for  our  proposed 
approach  which  we  named  ADCON  (for  Adaptive  Decentralized  CONtrol)  comes  from 
the  following  observations  about  the  current  status  of  control  theory. 

;  An  Important  aspect  of  centralized  control  has  been  the  study  of  systems 

*  with  unknown  or  uncertain  (time  varying,  random)  parameters.  The 

investigation  of  this  problem  led  to  an  extensive  literature  on  adaptive 
control  (also  called:  learning  or  self-organizing  systems).  The  natural 
progression  In  developing  centralized  controllers  was  from  the  non-adaptlve 
case  to  the  more  difficult  problems  addressed  by  adaptive  techniques. 

The  study  of  decentralized  control  seems  so  far  to  be  almost  exclusively 
devoted  to  non-adaptlve  techniques.  A  possible  explanation  of  this  state  of 
affairs  Is  the  fact  that  the  area  of  decentralized  control  of  completely  known 
i  systems  still  has  many  unresolved  Issues  and  some  basic  problems  are  yet  to  be 

answered.  Under  these  conditions,  there  seemed  to  be  little  Incentive  to 
;  tackle  the  more  complex  adaptive  case  which  deals  with  partially  known 

systems.  However,  this  line  of  thinking  Is  based  on  the  experience  gained  In 
centralized  control  and  It  may  be  Inapplicable  In  the  context  of  the 
decentralized  problem,  which  has  radically  different  characteristics.  In 
fact,  adaptive  techniques  have  a  central  role  In  decentralized  control,  which 
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Is  of  a  somewhat  different  nature  than  the  role  they  play  In  the  centralized 
problem. 

To  understand  the  Interrelation  between  adaptive  and  decentralized 
control,  we  have  to  re-examine  the  basic  Issues  underlying  the  need  for 
decentralized  control  strategies.  The  main  motivation  for  considering  such 
strategies  arises  In  the  context  of  complex,  large-scale  systems  where  a 
centralized  controller  usually  requires  excessive  computational  requirements 
and  excessive  Information  gathering  networks  to  make  such  a  controller 
feasible.  In  such  a  system.  It  Is  reasonable  to  assume  that  the  local 
controller  (l.e.,  the  controller  of  one  subsystem  In  the  large  system)  has 
only  partial  Information  about  the  rest  of  the  system.  Even  If  the  structure 
of  the  whole  system  (l.e.,  the  state  equations  of  all  subsystems  and  their 
Interactions)  can  be  made  available  to  each  local  controller,  the  sheer 
complexity  of  the  problem  often  limits  the  usefulness  of  this  Information.  In 
fact,  attempting  to  use  too  much  Information  may  be  one  of  the  principal 
stumbling  blocks  of  conventional  approaches  to  decentralized  control.  Most  of 
these  approaches  try  to  solve  the  (optimal)  centralized  problem,  and  then  to 
find  clever  ways  of  decentralizing  the  solution.  The  shortcomings  of  this 
technique  and  the  need  for  a  different  point  of  view  are  by  now  widely 
recognized. 

The  basic  Idea  underlying  our  approach  Is  to  assume  that  from  the 
subsystem's  point  of  view,  the  rest  of  the  system  Is  not  exactly  known.  Thus, 
the  subsystem  Is  aware  of  Its  own  structure,  but  It  has  only  an  approximate 
knowledge  of  the  rest  of  the  system,  for  example.  In  the  form  of  a  reduced 
order  model.  (Different  subsystems  will  use  different  models  of  the  "outside 
world".)  The  local  controller  Is  then  designed  on  the  basis  of  this  partial 
Information.  The  modeling  uncertainty  Inherent  In  this  procedure  makes  It 
necessary  to  consider  robust  or  adaptive  control  structures.  Note  that  the 
uncertainty  here  Is  due  to  the  complexity  of  the  system  rather  than  to  lack  of 
knowledge  or  to  random  effects,  which  are  the  traditional  sources  of 
uncertainty  In  centralized  control.  The  Idea  of  replacing  a  complex 
deterministic  problem  by  a  simple  stochastic  model  Is  by  no  means  new,  and  has 
been  used  In  a  variety  of  physical  problems  (e.g.,  statistical 
thermodynamics). 


The  use  of  reduced  order  models  and  partial  Information  greatly 
simplifies  the  design  and  Implementation  of  the  decentralized  controllers.  It 
raises,  however,  many  difficult  questions  regarding  the  conditions  under  which 
such  a  scheme  will  lead  to  satisfactory  system  behavior.  What  Is  needed  Is  a 
theory  for  the  control  of  Interconnected  subsystems  In  the  presence  of  model 
uncertainties.  In  an  earlier  report  [123  and  In  some  related  papers  we  made  a 
preliminary  study  of  some  of  these  Issues. 

An  even  more  difficult  set  of  questions  arises  with  regard  to  the 
operation  of  adaptive  controllers  In  the  presence  of  uncertainty.  Currently 
available  adaptive  control  algorithms  have  been  shown  to  experience  severe 
difficulties  In  the  presence  of  unmodeled  plant  dynamics.  We  were  able  to 
derive  conditions  which  guarantee  that  the  adaptive  controller  will  have 
specified  performance  despite  plant  uncertainty  and  unmodeled  dynamics.  These 
conditions  provide  guidelines  for  the  analysis  and  design  of  robust  adaptive 
controllers.  A  combination  of  results  from  robust  control  and  adaptive 
control  theory  was  used  to  prove  the  main  theorem.  The  main  theorem  was 
applied  to  a  number  of  well-known  adaptive  structures:  the  direct  adaptive 
controller,  an  adaptive  observer,  the  Indirect  adaptive  controller,  and  a 
general  form  of  the  model  reference  adaptive  controller  [4].  We  believe  that 
this  work  represents  a  significant  advance  In  the  field  of  adaptive  control. 

In  [13]  we  presented  an  Input-output  approach  for  analyzing  the  global 
stability  and  robustness  properties  of  adaptive  controllers  to  unmodeled 
dynamics.  The  concept  of  a  tuned  system  was  introduced,  l.e.,  the  control 
system  that  could  be  obtained  If  the  plant  were  known.  Comparing  the  adaptive 
system  with  the  tuned  system  results  In  the  development  of  a  generic  adaptive 
error  system.  Passivity  theory  was  used  to  derive  conditions  which  guarantee 
global  stability  of  the  error  system  associated  with  the  adaptive  controller, 
and  ensure  boundedness  of  the  adaptive  gains.  Specific  bounds  are  presented 
for  certain  significant  signals  In  the  control  systems.  Limitations  of  these 
global  results  are  discussed,  particularly  the  requirement  that  a  certain 
operator  be  strictly  positive  real  (SPR )  —  a  condition  that  Is  unlikely  to 
hold  due  to  unmodeled  dynamics. 

The  ADCON  concept  Involves  many  different  Issues,  as  can  be  seen  from  the 
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earlier  discussion  and  from  [4], [9], [12], [13].  So  far  we  have  addressed  the 
problem  of  designing  a  controller  for  a  single  subsystem,  when  the  rest  of  the 
system  Is  fixed.  This  represents  only  one  step  In  an  Iterative  procedure  In 
which  each  subsystem  performs  Its  own  controller  design.  We  have  done  some 
Investigation  extensions  of  the  theory  of  robust  control  and  adaptive  control 
to  the  case  of  Interconnected  subsystems,  In  which  local  controllers  are 
designed  sequentially  (Iteratively)  or  simultaneously.  A  number  of  different 
Information  structures  were  considered.  It  seems  that  by  providing  each 
subsystem  with  structural  Information  In  addition  to  an  aggregate  (reduced 
order)  model  of  the  rest  of  the  systems.  It  Is  possible  to  obtain  simpler 
design  schemes.  However,  no  conclusive  results  are  available  at  this  time. 

We  have  also  Investigated  the  application  of  lattice  structures  to  the 
adaptive  control  problem.  Our  work  In  this  area  seemed  to  have  generated  a 
considerable  amount  of  Interest  (cf.  [R1]-[R6] ) .  This  class  of  algorithms  Is 
especially  well  suited  for  large  scale  problems  of  the  type  considered  In  this 
project. 

In  the  next  section  we  list  the  publications  prepared  under  this 
contract.  The  key  papers  are  enclosed  In  the  appedlces. 
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1.  INTRODUCTION 


1.1  Background 

The  analysis  and  design  of  adaptive  control  systems  has  been  the  subject 
of  extensive  research  in  the  past  two  decades  [1]-[10].  Adaptive  techniques 
provide  a  way  of  handling  plant  uncertainty  by  adjusting  the  controller 
parameters  on-line  to  optimize  system  performance.  An  alternative  method  for 
handling  uncertainty  is  to  use  a  fixed  structure  controller  designed  to 
provide  acceptable  performance  for  a  specified  range  of  plant  behavior.  In 
principle,  adaptive  controllers  can  provide  Improved  performance  compared  to 
fixed  robust  controllers,  since  they  are  tuned  to  the  uncertain  plant. 

However,  adaptive  controllers  sometimes  exhibit  undesirable  behavior  during 
the  tuning  or  adaptation  process.  For  example,  unmodeled  dynamics  can  cause  a 
rapid  deterioration  In  performance  and  even  instability  [11], [12].  This 
problem  Is  not  resolved  by  Increasing  the  order  or  complexity  of  the  model. 
Since  the  model  of  any  dynamic  system,  by  definition.  Is  not  the  actual 
system,  it  can  therefore  be  argued  that  unmodeled  dynamics  are  always  oresent, 
ad  infinitum. 

The  main  reason  for  these  difficulties  with  adaptive  controllers  seems  to 
be  that  robustness  to  unmodeled  dynamics  was  not  considered  as  a  design 
criterion  in  the  development  of  the  adaptive  control  algorithm.  The  design 
objective  Is  global  stability  of  the  closed-loop  system,  e.g.,  [7],  [9]  and 
various  assumptions  on  the  structure  of  the  plant  are  required  to  achieve  that 
objective.  In  particular.  It  is  necessary  to  assume  that  the  plant  Is  linear 
and  time  Invariant  (LTI),  that  the  relative  degree  of  the  transfer  function  Is 
known  as  well  as  the  sign  of  the  high  freauency  gain.  Such  requirements  are 
not  practical  since  real  plants  are  often  nonlinear  and  time-varying  and  can 
be  accurately  represented  only  by  high  order  (sometimes  infinite  order  [13]) 
complicated  models. 

The  need  for  robustness  to  plant  uncertainty  is  not  unique  to  adaptive 
control.  The  problem  of  robustness  Is  ubiquitous  in  control  theory  and  has 
been  studied  in  the  context  of  fixed  (nonadaptlve)  control  [14]— [17] .  These 
studies  rely  on  the  Input/output  properties  of  systems,  e.g.,  [18], [19].  The 
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predominant  reason  to  examine  robustness  Issues  In  this  way  Is  that  the 
characterl tl cs  of  unmodeled  dynamics,  such  as  uncertain  model  order,  are 
easily  represented.  Lyapunov  theory,  on  the  other  hand.  Is  not  well  suited 
for  this  type  of  uncertainty.  Typically,  plant  uncertainty  Is  characterl zed 
by  assuming  that  the  plant  belongs  to  a  well  defined  set.  For  example,  a  set 
description  of  an  uncertain  LTI  plant  Is  to  define  a  "ball"  In  the  frequency 
domain.  The  center  of  the  ball  Is  the  nominal  plant  model,  and  the  radius 
defines  the  model  error.  This  set  model  description  is  one  type  of  a  more 
general  set  description,  referred  to  as  a  conic- sector  [15].  The  uncertainty 
In  the  plant  Induces  an  uncertainty  in  the  input/outout  map  of  the  closed-loop 
system  which  can,  again  be  characterl zed  by  a  conic  sector.  Performance 
requirements  for  the  control  system  can  be  translated  into  statements  on  the 
conic  sector  which  bounds  the  closed-loop  systems,  making  It  possible  to  check 
whether  a  given  design  meets  specifications,  and  providing  guidelines  for 
robust  controller  design. 
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In  this  paper  we  use  the  Input/output  approach  to  analyze  the  global 
stability  and  robustness  properties  of  contl nuous-tlme  adaptive  controllers 
with  respect  to  unmodeled  dynamics  (although  we  consider  only  contl nuous-tlme 
algorithms,  the  Input-output  formalism  can  be  readily  extended  to  the 
discrete-time  case).  By  global  we  mean  that  no  specific  magnitude  constraint 
(other  than  boundedness)  Is  placed  on  any  of  the  external  Inputs  or  Initial 
condltons.  We  develop  an  adaptive  error  system  of  a  general  form,  by 
comparing  the  actual  adaptive  system  with  a  tuned  system,  l.e.,  the  control 
system  that  could  be  obtained  If  the  plant  were  known.  This  error  system  Is 
similar  to  the  type  used  In  [7], [8]  where  the  tuned  system  error  output  is 
zero,  due  to  the  assumption  of  perfect  modeling.  By  relaxing  this  assumption 
we  show  that  the  non-zero  outputs  of  the  error  system  are  the  Inputs  to  a 
nonlinear  feedback  error  system  consisting  of  the  adaptive  algorithm  and  two 
feedback  (Interconnection)  operators, denoted  by  Hey  and  Hzy  . 


J  ' 


An  Important  consequence  of  this  structure  Is  that  the  existence  of 
solutions  (e.g.,  tuned  system  performance)  Is  separated  from  the  stabllty 
analysis  (e.g.,  stability  of  the  nonlinear  error  system).  In  general,  the 
adaptation  law  Is  passive;  consequently.  If  Hev  Is  strictly  positive  real 
(SPR),  then  application  of  passivity  theory  [19]-[21],  provides  global 


Instability  of  the  map  from  the  tuned  system  output  to  the  actual  adaptive 
system  output,  even  though  the  adaptive  parameters  may  grow  beyond  all 
bounds.  We  provide  other  conditions  (e.g.,  Hzv  stable)  to  Insure  the 
l  boundedness  of  the  adaptive  gains.  Similar  results  are  developed  to 
Insure  L  -stability  of  the  error  system  by  using  an  exponentially  weighted 

m 

passivity  theory  [19].  These  results  are  summarized  In  Theorems  1A  and  IB. 

As  a  by  product  of  the  Input/output  view  we  also  obtain  specific  bounds 
on  the  I2  and  norms  of  significant  signals  In  the  adaptive  system.  The 
results  are  summarized  In  Corollary  1. 

The  results  In  Theorem  1  and  Corollary  1  are  not  essentially  new  (see 
e.g.,  [7], [8]),  although  they  do  provide  some  extentions  to  previous 
results.  The  main  contribution,  however.  Is  the  fact  that  all  the  results  can 
be  obtained  from  a  generic  error  system  and  from  the  application  of  nonlinear 
stabllty  theorems  based  on  Input-output  properties.  As  a  consequence  of  this 
approach.  It  is  to  be  expected  that  conditions  for  robustness  will  arise  In  a 
natural  way.  Such  robustness  results  are  obtained,  but  unfortunately, they 
have  a  limited  practical  use.  The  main  limitation  Is  that  the  global  theory 
(Theorem  1)  requires  that  Hey  e  SPR  ,  which  In  turn  places  an  upper  bound  on 
the  size  of  the  unmodeled  dynamics  In  the  plant.  The  details  are  contained  In 
Lemmas  4.1  and  5.2.  This  bound  Is  quite  restrictive  and  Is  easily  violated  by 
even  the  most  benign  model  errors,  thus,  verifying  the  results  obtained  In 
[11],  [12].  To  over  come  this  limitation,  we  construct  an  SPR  compensator, 
based  on  the  scheme  proposed  In  [22]  In  the  context  of  robust  (non-adaptlve) 
control.  Although  in  the  adaptive  case  the  supporing  arguments  are  heuristic, 
an  example  simulation  shows  a  positive  result. 

The  Input/output  analysis  presented  here  provides  a  generic  framework 
within  which  It  Is  possible  to  analyze  the  robustness  of  adaptive  robust 
controllers.  We  believe  that  this  framework  can  be  used  to  develop  practical 
adaptive  control  algorithms  that  can  be  more  readily  applied  to  real  systems, 
than  the  class  of  algorithms  currently  In  use. 

Since  this  paper  merges  Ideas  from  several  areas.  It  is  necessary  to 
Introduce  a  number  of  definitions  and  concepts. 


Since  this  paper  merges  Ideas  from  several  areas.  It  Is  necessary  to 
Introduce  a  number  of  definitions  and  concepts. 


2.  SOME  PRELIMINARIES 


2.1  Hotatl on 

The  Input/output  formulation  of  multivariable  systems  Is  the  principal 
view  taken  throughout  this  paper  and  the  notation  and  terminology  used  is 
standard  (see  e.g.  [18] ,[19]).  The  Input  and  output  signals  are  assumed  to  be 
Imbedded  in  either  the  normed  function  space 

ij;  ■  (X  :  CO,.)  .  Rn  |  ||,||P  <  -(  (2.U) 

or  Its  extentlon 

Lpe  «  (x  :  [0,T]  ♦  Rn|  ||x||Tp  <  -,  T  <  -}  (2.1b) 

The  respective  norms  ||*|lp  and  j|.||Tp  are  defined  as  follows: 

1 1*1  lp  •  Hw  11*1  lTp  (2.2a) 

with 

/(/  |x(t)|pdt)l/p  ,  P  c  Cl,-) 

l|x|lTp  •  |  °  (2.2b) 

(  _  Ix(t)|,  p  ■  - 

V  t*[0,T] 

where  |.|  Is  the  Euclidean  norm  on  Rn.  Hence,  Is  an  Inner  product 
space,  with  inner  product  <x,y>T  of  elements  x,  y  c  lge  defined  by 

T 

<x,y>_  ■  /  x(t)'y(t)dt  (2.3) 

o 

and  so  jjx||j2  •  (<x,x>T)1^2  •  If  T  ♦  •  then  an  Inner-product  space 

with  Inner  product  <x,y>  »  11m<x,y>  . 

T-m.  1 


2.2  Stablllt 


Systems  considered  In  this  paper  are  described  by  Input/output  equations 

of  the  form  y  ■  Gu  where  G:L™  ♦  l"  Is  a  causal  map  from  u  Into  y,  also 

pe  pe 

denoted  u  ♦  y  •  The  system  G  Is  said  to  be  Lp-stable  (or  simply  stable)  If  G 

maps  u  e  L™  Into  y  e  Ljj  and  If  there  exists  finite  constants  k  and  b  such 

that  | jGu] |Tp  <  k  ||u|lTp  ♦  b  ,  for  all  T  >  0  and  all  ueLpe  •  The  smallest 

k  that  can  be  found  Is  referred  to  as  the  Lp-ga1n  (or  simply  pain)  of  G, 

denoted  y„(G)  • 

P 

Because  we  often  encounter  LTt  systems  It  Is  convenient  to  introduce  the 

following  notation.  Let  R(s)  and  RQ(s)  denote  the  proper  and  strictly  proper 

rational  functions,  respectively.  Let  S  and  S  denote  functions  In  R(s)  and 

o 

R^( s)  ,  respectively,  whose  poles  all  have  negative  real  parts.  Thus, 

S  and  Sq  are  the  stable,  lumped,  LTI  systems.  Denote  multivariable  systems 

with  transfer  function  matrices,  by  R(s)nxm,  Snxm  .  etc.  For  example, 
nxm 

G  e  SQ  means  that  all  elements  of  G  belong  to  SQ  ,  and  so  on. 

If  G  e  Snwn  then  the  following  Lp-ga1ns  are  obtained. 


Y, (G)  <  y  (G)  »  f  7[G(t)]dt 

1  *  0 

(2.4) 

Y,(G)  *  sup  TtG(ju)] 

(2.5) 

where  7(A)  denotes  the  maximum  singular  value  of  the  matrix  A,  defined  as  the 
positive  square  root  of  the  maximum  eigenvalue  of  A* A,  where  *  Is  the 
conjugate  transpose  of  A.  In  (2.4),  (2.5)  G  is  the  operator,  G(j«)  the 
transfer  function  matrix,  and  G(t)  is  the  Impulse  response  matrix. 

2.3  Passivity 

The  following  definitions  follow  those  In  [19], [21].  Let 
G;L™e  ♦  *nd  let  u,  p  be  constants  with  u  >  0  .  Then,  Vue  L™e  : 


6  ts  gasslve  If, 

<  u,  G  u>j  >  p  (2.6) 

G  Is  Input  strictly  passive  If, 

<  u,  Gu  >T  >  p  ♦  u|UI-f2  (2.7a) 

G  Is  output  strictly  passive  if, 

<  u,  Gu  >j  >  p  *  tiiGuij^  (2.7b) 

(u  and  p  are  not  the  same  throughout).  When  G  e  Sraxm  satisfies  (2.7),  G  Is 
said  to  be  strictly  positive  real  (SPR),  denoted  G  e  SPR™  .  Because  SPR 
systems  play  a  crucial  role  In  the  proof  of  stability  of  adaptive  systems,  we 
Introduce  the  following  subsets: 

SPR™  •  (G  c  S™*™!^  [G(  ju)  ♦  G( -  ju ) '  ]  -  pi)  >  0,  VucR}  (2-8«) 

SPR™  '  (G  s  [G(ju)  ♦  G(-ju)1]  -  p  G( -ju) 'G( ju))  >  0,  VueR}  (2.8b) 

where  x_( A)  denotes  the  smallest  eigenvalue  of  A.  Thus,  whenever  g  e  S™x™  . 
conditions  (2.7)  can  be  tested  In  the  frequency  domain.  Moreover,  SPR™  and 
SPR™  ,  respectively,  separate  the  strictly  proper  SPR  functions  from  the 
proper,  but  not  strictly  proper,  SPR  functions.  In  the  scalar  case,  the 
frequency  domain  conditions  simplify  because  aTGUu)  +  G(-ju)']* 

2  ReCG(ju)]. 

Certain  unstable  systems  In  R(s)mxm  can  be  passive  by  virtue  of  (2.6). 

In  particular,  GeRfs)™5*™  Is  passive  If  G(s)  Is  positive  real.  The  transfer 
function  matrix  G(s)  Is  positive  real  If:  (1)  It  has  no  poles  In  Re(s)  >  0, 
(11)  poles  on  the  ju  axis  are  simple  with  a  non-negative  residue,  and  (111) 
for  any  u  c  R  not  a  pole  of  G(jw)  +  G(-Ju)'  >  . 


Si 
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2.4  Model  Error 


The  cornerstone  of  robust  control  design  Is  a  quantifiable  bound  on  the 
error  between  the  model  used  for  control  design  and  the  actual  plant  to  be 
controlled.  In  the  adaptive  control  case  considered  here  the  model  Is  a 
parametric  model,  where  the  parameters  are  not  known  exactly.  The  structure 
of  the  parametric  model  can  be  obtained  analytically  from  physical  laws,  but 
this  Invariably  results  in  a  complicated  model.  Often  a  simple  structure  Is 
selected  because  It  Is  more  convenient  for  analysis  and  synthesis. 


Let  P  denote  the  plant  to  be  controlled.  In  the  broadest  sense  P  Is  a 
relation  In  L?_  x  l?  ,  l.e.,  the  set  of  all  possible  ordered  pairs 


*  Li*  x  Lit 

(u,y)eL™e  x  l"e  of  Inputs  u  e  L™e  and  outputs  yeL”e  that  could  be  generated 
by  the  plant  [18].  The  uncertainty  in  the  plant  is  denoted  by  (u,y)  e  P  . 


Let  P  tL91  ♦  Ln  denote  a  parametric  model  of  the  plant  P  with 
a  pe  .  pe 

parameters  a  e  R*  .  The  parameters  can  be  selected  so  as  to  minimize  any 
discrepancies  between  the  model  and  the  plant,  l.e.. 

Inf iy-PauiTp  -  iy-P*uiTp  (2.9) 

aeR  * 


We  will  refer  to  ^.^r*  as  the  tuned  model  parameters  and  to  P  ■  p^.  as  the 

- -  a# 

tuned  parametric  model  of  the  plant.  In  general,  P*  Is  dependent  on  the 
Input/output  sequence. 

Most  of  the  previous  work  on  adaptive  control  deals  with  the  case  where 

for  every  (u,y)  e  P  there  exists  a  tuned  parametric  model  P*,  such  that 

P**P.  In  this  paper  we  consider  the  presence  of  unmodeled  dynamics,  thus, 

the  uncertain  plant  P  cannot  be  perfectly  modeled  by  any  parametric  model 

P  .  Since  we  will  deal  exclusively  with  LTI  plants  P  e  R(s)nxm  ,  It  Is 
a 

convenient  to  describe  this  model  error  In  the  frequency-domain.  Let 
Bs(r)  denote  a  "ball"  In  S  of  radius  r,  defined  by 


B$(r)  (6  c  Snxm|  ^G(ju)]  <  r<«),  «  c  R} 


(2.10) 


Let  the  plant  to  be  controlled  be  described  by 


P  «  (1  +  a)P* 


(2.11a) 


where  P  e  R(s)nxm  Is  the  plant,  p*  e  R(s)nxm  Is  the  tuned  parametric  model, 
and  a  e  Snxn  denotes  the  unmodeled  dynamics.  Further,  the  only  knowledge 
available  about  a  is  that  It  Is  bounded  such  that 

a  c  B$U)  (2.11b) 

where  6(u)  Is  known  for  all  frequencies.  In  other  words,  while  the  operator 
a  Is  not  precisely  known,  we  do  know  a  bound  on  Its  effect.  This  model 
description  (2.2)  Is  used  throughout  the  paper  to  precisely  define  the  plant 
to  be  controlled  in  an  adaptive  system.  Following  Doyle  and  Stein  [16]  we 
will  refer  to  (2.11b)  as  an  unstructred  uncertainty.  Mote  that  although  a  Is 
stable,  P  and  P*  need  not  be  stable.  Hence,  the  parametric  model  Is 
Implicitly  required  to  capture  all  unstable  poles  of  the  plant.  Although  this 
is  not  severly  restrictive  -  at  least  on  practical  grounds  -  nonetheless.  It 
can  be  eliminated  by  deflnng  model  error  as  (stable)  deviations  In  (stable) 
coprime  factors  of  the  plant  [23].  As  the  subsequent  analysis  is  not 
substantially  effected  by  this  choice,  we  will  remain  with  (2.11)  for  purposes 
of  Illustration. 

2.5  Persistent  Excitation 

From  [31],  a  regulated  function  F(.)  ■  R+  ♦  Rnxni  Is  persistently 
exciting,  denoted  F  e  PE  ,  If  there  exists  finite  positive  constants 
a^,  ag,  and  such  that 

*+a3 

a-  l  >  f  F(t)F(t)'dt  >  a,  I  ,  »  S  e  R,  (2.12) 

z  n  j  in  ♦ 

The  usefulness  of  a  persistently  exciting  signal  Is  In  establishing  the 
exponential  stability  of  the  following  differential  equation  which  arises  In 
many  adaptive  and  Identification  schemes,  l.e.. 


x  «  -BFHF’x  ♦  m  , 


x{0)  e  R 


(2.13) 


It  is  shown  In  [31]  that  If  Be  Rnxffl,  B  «  B*  >  0  ,  He  SPRj  or  SPRJ,  and 
F  e  PE  ,  then  (w,  x(0)f-*x  Is  exponentially  stable,  l.e.,  3  m,  x  >  0  such 

that 


I x( t) |  <  me'xt  (x(0) |  ♦  /  me“x{t_r)  |w(T)|dT  • 

0 

We  will  utilize  this  latter  result  In  section  IV  In  our  proof  of  stablity  of 
the  adaptive  system. 


% 

% 
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3.  ADAPTIVE  ERROR  MODEL 

In  this  section  we  develop  a  generic  adaptive  error  model  which  will  be 
used  in  the  subsequent  analysis.  This  requires  defining  the  notions  of  robust 
control  and  tuned  control. 

Robust  and  Tuned  Control 


Consider,  for  example,  the  model  reference  adaptive  control  (MRAC) 
depicted  In  Figure  3.1,  consisting  of  the  uncertain  plant  P,  a  reference  model 
Hr,  and  an  adaptive  controller  C(a)  ,  where  e  is  the  adaptive  gain  vector,  r 
is  a  reference  Input,  d  Is  a  disturbance  process,  and  n  Is  sensor  noise. 

Oenote  by  H(e)  the  closed-loop  system  relating  the  external  inputs  w  *  (r' , 
d',  n')'  to  the  output  error  e,  as  depicted  In  Figure  3.2..  Also,  let  w  e  W 
denote  the  admlssable  class  of  Input  signals. 

The  objective  of  the  adaptive  controller  Is  twofold:  (1)  adjust  a  to  a 
constant  a*  e  Rk  such  that  H(e*)  has  deslreable  properties;  and  (2)  during 
adaptation,  as  a  is  adjusted,  the  error  Is  well  behaved.  In  the  usual 
formulations  [7]  only  (1)  Is  considered  and  further  It  Is  assumed  that  there 
exists  a  matched  gain,  denoted  by  3"  e  R  .  such  that 

H(TT)  ■  0  (3.1) 

The  presence  of  uncertain  unmodeled  dynamics  In  the  plant  eliminate  the  chance 
of  satisfying  the  matching  condition.  Thus,  It  is  more  appropriate  to  define 
a  tuned  gain,  denoted  by  a*  e  R*1  ,  corresponding  to  each  (u,y,w)  e  P  x  W  , 
such  that 

H(a«)w  <  H(e)w  ,  v  a  c  Rk  (3.2) 

The  error  signal  e*  :■  H(a*)w  Is  referred  to  as  the  tuned  error.  Note  that 
each  (u,y,w)  e  P  x  W  engenders  a  possibly  different  a*  .  Also,  It  is 
Important  to  distinguish  the  tuned  gain  ,  from  the  robust  gain  a.  c  Rk  . 

w  0 

where 

sup  H(e  )w  <  sup  H(a)w  ,  V  a  e  R*  (3.3) 

P  x  M  °  P  x  W 
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The  error  signal  e  :«  H(e  )w  is  referred  to  as  the  robust  error.  It  follows 
o  o  ■ 

from  these  definitions  that  the  tuned  error  is  always  smaller  in  norm  than  the 

robust  error,  thus  V  w  e  W  , 

e*  *  H(e*)w  <  eQ  *  H(8q)w  ,  (3.4) 

The  tuned  controller  is,  unfortunately,  unrealizable  since  it  requires  prior 
knowledge  of  the  actual  system  H(a)  (or  equivalently,  the  plant  P)  and  the 
input  w.  A  practical  adaptive  controller  is  likely  to  have  a  larger  error 
norm. 

Structure  of  the  Adaptive  Control 

In  summary,  we  consider  the  multivariable  adaptive  system,  shown  In 
Figure  3.2,  and  described  by 


e  *  H(a)w  .  (3.5) 

where  e(t)  e  Rm  is  the  error  signal  to  be  controlled,  w(t)  e  Rq  Is  the 

*  if 

external  input  restricted  to  some  set  W,  and  e(t)  e  R  Is  the  adaptive 
gain.  The  class  of  adaptive  controllers  considered  here  are  such  that  the 
adaptive  gains  multiply  elements  of  Internal  signals  z(t)  e  R  ,  referred  to 
as  the  regressor,  to  produce  the  adaptive  control  signals, 

ft  ■  81  Zj  ,  1  e  Cl.m]  (3.6) 

where  And  z^  are  k^-dlmenslonal  subsets  of  the  elements  In  a  and  z, 
respectively.  Thus, 

k  «  i  k.  (3.7) 

1«1  1 

Oef 1 ne  the  adaptive  gain  error, 

e(t)  :«  a(t)  -  a*  (3.8) 

where  9*  c  R  Is  the  tuned  gain  (3.4).  Also,  define  the  adaptive  control 
error  signals. 


f  f 


Vj  :«  9}  z^  ,  1  «  1,  ....  m 


(3.9) 


An  equivalent  expression  Is, 


v  =  Z'0 


where  the  time-varying  matrix  Z  Is  defined  by 


(3.10a) 


Z  *  block  d1ag(zlt  z2,  .  .  .,  z*) 


(3.10b) 


To  describe  the  relations  among  the  signals  e,  z,  v,  and  w  we  Introduce 
the  Interconnection  system  :  (w,v  )  ♦  (e,z)  ,  as  shown  In  Figure  3.3.  In 

particular,  let  Hj  e  R(s)^m+lc^m+q^  t  and  wnere  Hj  Is  defined  by. 


:«  H, 


Hew  \  /  “Hev  w 


(3.11) 


In  effect,  this  structure  serves  to  Isolate  the  adaptive  control  error  v,  from 
the  rest  of  the  system.  When  the  adpatlve  control  Is  tuned,  9*0  and  v  «  0; 
consequently,  the  tuned  error  signal  (3.4)  Is, 


e*  H<8jw  *  H  w 
*  *  ew 


We  can  also  define  a  tuned  regressor  signal. 


z*  :■  H  w 
*  zw 


(3.12) 


(3.13) 


In  general,  all  the  subsystems  In  Hj  are  dependent  on  the  tuned  gains  a* 


The  Interconnection  system  can  also  be  written  as. 


.■•V*  .•  .*•  /*  ►  - .  • .  *  *  *  •’*  '  *  v  ■.*  >  >*%***.'•  v  v  *«*** 


e  •  e*  -  H^v  (3.14a) 

2  *  2*  -  Hzv  v  (3.14b) 

with  v  given  by  (3.10).  To  complete  the  error  model  requires  describing  the 
adaptatlve  algorithm.  I.e.,  the  means  by  which  e(t)  Is  qenerated.  Me  will 
consider  two  typical  algorithms.  A  constant  gain  (gradient)  algorithm  [7]: 

• 

9  *  r  Z  6  (3.15) 


where  r  e  Rkx,c,  r  *  r'  >0  ,  and  a  similar  but  nonlinear  gain  algorithm: 

o  *  r(Ze  -  p(e)e)  (3.16a) 


where  p  :  Rk  ♦  R+  Is  a  retardation  function,  whose  purpose  is  to  prevent 
e  from  growing  too  quickly  in  certain  situations.  Although  many  functions 
will  suffice  we  will  select  the  one  proposed  In  [24],  namely: 


(m/c  -  l)z,  i0i  >  c  :«  maxi 9*i 


>  ( 9  ) 


,  iei  <  c 


(3.16b) 


The  complete  adaptive  error  system.  Is  shown  In  Figure  3.4.  Note  that 
the  error  system  Is  composed  of  two  subsystems:  a  linear  subsystem  el  and  a 
non-linear  subsystem  z  . 


ijy 
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4.  CONDITIONS  FOR  GLOBAL  STABILITY 


The  theorems  stated  below  give  conditions  for  which  the  adaptive  error 
system  (Fig.  3.4)  Is  guaranteed  to  have  certain  stability  and  performance 
properties.  Proofs  are  given  In  Appendix  A.  Heurlstlcally,  however,  the 
basis  for  the  proofs  Is  application  of  the  Passivity  Theorem  ([193.  pg.  182). 
It  turns  out  that  the  map  e  ♦  v  Is  passive.  Thus,  If  Hev  Is  SPRra  ,  then 
the  map  ft,  ♦  (e,v)  Is  L2*stab1e  even  though  z  and/or  e  can  grow  without 
bounds.  Further  restrictions,  provided  below,  cause  9  and  z  to  be  bounded. 
(We  use  the  notation  Hx  ♦  0  (exp.)“  to  mean  that  x(t)  ♦  0  (exponentially)  as 
t  ♦  •  . ) 

Theorem  A:  Global  Stability 


For  the  adaptive  error  system  shown  in  Figure  3.4,  assume  that: 

(Al)  The  system  Is  well -posed  In  the  sense  that  all 
Inputs  w  e  w  produce  signals  e.v.z,  e  ,  and 

8  In  L  . 

•e 

(4.1a) 


( A2) 

ckxm 

Hzv  e  So 

(4.1b) 

(A3) 

Hey  c  SPR™ 

(4.1c) 

Under  these  conditions: 

(1) 

If  (e*.  e*)  e  L*  Hl"  H4  e,  -*-0)  and  (z*,  z,)  c  L* 
algorithm  (3.15)  or  (3.16): 

then  with 

( 1-a)  (e.e)  t  lJ,  e  t  Lg  n L*  ,  and  e  -►0. 

(4. a) 

(1-b)  e  6  L?nL!*  *  e  LH»  and  -**0- 

(4.2b) 

(1-c)  v  c  L™n  lJ  ,  y  t  L™  ,  and  v  ♦  0  . 

(4.2c) 
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t 


(1-d)  (z.i)  e  t*  ,  (z-z*,  l-ij  e  l\ni*  .  and  z-z*-*0  exp. 

(4. 2d) 

(1-e)  If,  In  addition,  e#  ■  0  (matched)  and  z*  e  PE  then 
(e,  9,  e-e*,  v,  z-z*)  —►0  exp. 

(4.2e) 

If  (e*,  e*)cL™  and  (z*,  z*)  e  i*  ,  then  with  algorithm  (3.15): 


i 


(1i-a)  ZclJ  (4.3) 

(11-b)  With  the  addition  of  either  algorithm  (3.16)  or  z  e  PE  It  follows 
that  the  elements  of  a,  8,  e,  e,  v,  v,  and  z  are  In  . 

(4.4) 

Theorem  IB:  Global  Stability 

Replace  (A3)  In  Theorem  1  by 


(A3)‘  HevcSPRS 

(4.5) 

If  (e*.  e*)  e  ijni™  («>  e*-*0)  ,  and  (z*. 

IJ  e  L*  then  with 

algorithm  (3.15)  or  (3.16) 

n-a)  (a,  a)e  l*  ,  e  e  l^nL*  ,  e-m-0 

(4.6a) 

(1-b)  edjnu;  ,  ee  l;  ,  e  -  e*-*0 

(4.6b) 

(1-c)  ( v, v)  c 

(4.6c) 

H-d)  (z,z)  c  L *  ,  (z-z*,  z-z*)  e  ijjH  i*  , 
and  z-z^—^Q  . 

<4.6d) 

(1-e)  If,  In  addition,  e*  ■  0  (matched)  and  z#e  PE  , 

then  (a,  v)— ►O  exp.  (4.6e) 


•f.  t, 
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(11)  If  (e*,  fe*)  e  l"  and  {2*,  i*)  «  ,  then  with  algorithm  (3.15): 


(11-a)  z  c  Lk 


(4.7d) 


(11-b)  With  the  addition  of  either  zcPE  or  algorithm  (3.16),  the 
elements  of  8,  e,  e,  e,  v,  v  ,  and  2  are  In  L^. 


(4.7b) 


Corollary  1:  Performance  Bounds 


Suppose  z*  and  e*  satisfy  the  conditions  in  (1)  of  Theroems  1A  or  IB. 

(I)  Let  Hfiv  e  SPR*  ,  i.e.,  3  u.  y  >  0  such  that  VweR  , 

7CHev(jw)]  <  y  and  £[Hev(ju)  ♦  H^t-jw)']  >  w  Im  (4.8a) 

Then,  bounds  on  iei„  and  i9i  can  be  obtained  from: 

z  • 

ie-e«i2  <  Jjj-  [ie*i2  ♦  (ie*i2  +  2U  e(0)*  r_1o(0))1/2]  (4.8b) 

le'r^ei,,  <  e(0) *  r-19<0)  +  2iei2ie-e*i2/Y  (4.8c) 

(II)  Let  Hev  e  SPRq  .  1.e.,3»»  0.  k  >  0  such  that  V  u  e  R  . 

^(HevU«)  ♦Hev(-j«)*]  >w  Hev(-jtt)*  Hev(ju)  (4.9a) 

^Gey(4«)  ♦  Gev(-ju»)*]  >  k  lm  (4.9b) 


Gey(s)  :«  (1  ♦  qs)  H@v(s) 


(4.9c) 


Then,  bounds  on  ie.i  and  tei  can  be  obtained  from: 

z  • 

««i2  «  ^ie*+qe*i  ♦  (ie*+qe*i2  ♦  2k2U9(0) *r"1e(0))1/2] 
le'r^ei^  <  9(0)'r’le(0)  +  i  ie«+qe*i2iei2 


(4.9d) 


(4.9c) 


( 


I 


■ 

I 


I. 


Discussion 


(1)  Theorems  1A  end  18  give  conditions  under  which  the  adaptive  error 
system  Is  globally  stable.  Essentially,  conditions  are  Imposed  on  the 
Interconnection  subsystems  In  Hj  .  In  particular,  c  SPRm  and 

H  6  C  are  ^rect  retirements,  whereas  the  restrictions  on  the  tuned 
signals  e„  and  z*  ,  Indirectly  impose  requirements  on  and  .  These 
latter  requirements  are  dependent  on  knowledge  about  w  e  W  .  For  example.  If 
w  is  a  constant,  then  the  assumption  that  e*  ♦  0  (Theorem  1A-1)  requires 
that  the  tuned  feedback  system  is  a  Type-I  robust  servomechanism,  l.e.,  the 
transfer  junction  H^tO)  *  0  for  all  (u,y)  e  P  . 

(2)  Corollary  1  gives  explicit  bounds  on  signals  in  the  error  system. 

These  bounds  can  be  used  to  evaluate  the  adaptive  system  design.  Moreover, 

the  bounds  allow  a  coarse  determination  as  to  the  efficacy  of  adaptive  control 

vs.  robust  control.  By  comparing,  for  example,  the  adaptive  error  lei^  from 

(4.8)  with  the  robust  error  ie  i_  from  (1.5),  It  Is  possible  to  obtain  a 

o  z 

quantifiable  measure  of  performance  degradation  during  adaptation. 

(3)  Although  Theorems  1A  and  18  are  essentially  the  same,  there  are 
slight  difference  worth  noting.  These  differences  arise  because  In  1A, 

HeveSPR™—»  Hey(s)  Is  proper  but  not  strictly  proper,  whereas  In  IB, 
HeveSPR0"*  Hev^  is  str1c^y  proper.  Thus,  comparing  part  (1)  In  1A  and 
IB,  we  see  that  In  18,  v,  v  e  l”  whereas  In  1A,  v  is  addl tonally  In 
and  v-m»0  . 

(4)  The  use  of  persistent  excitation  or  gain  retardation  Is  seen  in  part 
(11)  of  theorems  1A  and  IB  to  provide  the  means  to  guaranty  bounded  signals. 
Other  schemes  based  on  signal  normalizations  or  dead-zones  can  provide  similar 
results,  e.g.  [32], [33].  The  effect  of  these  conditions  Is  to  provide  an 

L^-stablllty  which  Is  not  present  otherwise.  The  persistent  excitation 
condition  actually  supplies  exponential  stability,  which  Is  stronger  than 

l^-stablllty,  as  provided,  for  example,  by  the  gain  retardation  (see  proof  In 
Appendix  A). 

(5)  The  persistent  excitation  requirements  In  parts  (1)  and  parts  (11) 
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are  different.  In  parts  (1),  z#ePE  .  whereas  In  parts  (11),  tePE  •  The 
different  assumptions  arise  because  In  parts  (1)  we  enforce  the  matched 
condition  e^*0  .  Hence,  z*ePE  ■>  zePE  .  This  follows  from  (1-d)  where 
z  -  z#  ♦  0  expoentlally.  Also,  with  e*  ■  0  ,  a  bounded  disturbance  added  to 
the  reference  can  cause  z  c  PE  without  forcing,  e*  e  l  .  In  parts  (11), 
which  Is  more  realistic,  we  disallow  the  matched  condition,  and  hence, 

e  L^.  Thus,  z  e  PE  Is  the  weakest  assumption  to  make.  However,  since  z 
is  inside  the  adaptive  loop,  it  is  very  different  to  guarantee  z  e  PE  by 
Injecting  external  signals.  Note  also  (in  both  parts(ii))  that  without 
retardation  or  PE  it  is  possible  for  the  regressor  to  remain  bounded  even 
though  the  adaptive  parameters  may  grow  unbounded.  Similar  results  have  been 
reported  elsewhere,  e.g.  [24]. 

Robustness  to  Unmodeled  Dynamics 

Since  the  theorems  Impose  requrlements  on  the  Input/output  properties  of 

the  Interconnection  system.  It  follows  that  the  effect  of  model  error  on  these 

properties  determines  the  stability  robustness  of  the  adaptive  system.  For 

example,  both  theorems  reoulre  that  H_  e  SPRm  .  Suppose,  however,  that 

ev 

H  has  the  form, 
ev 

H  ■  (I  +  fl  )H  (4.10) 

ev  ev  ev 

where  H  Is  the  projection  onto  H  of  the  plant  uncertalny  operator  a  ; 

_  ev  ev 

and  Hgv  is  the  nominal  transfer  function  when  there  Is  no  uncertainty,  l.e., 

when  a  ■  0  .  Thus,  iT  Is  a  function  of  the  tuned  parametric  model  P*  and 

ev 

the  tuned  controller  gains  8*  .  (See  Section  Y  for  more  specific  formulae, 
e.g.  (5.S) .) 

Conditions  to  Insure  that  H  e  SPR?  despite  uncertainty  In  H„v  is 

ev  + 

provided  by  the  following: 

Lemma  4.1:  Let  Hev  be  given  by  (4.3).  Then  Hfty  e  SPR™  If  the  following 
conditions  hold: 

(1)  H#y  c  SPRj  (4.11a) 


*jr 


V. 

V 

•a 


Sr; 


v* 
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(4.11b) 


(11)  Hevc  8s(k)  where  V  «  e  R  , 

Mu.)  <  7  x[Hey(jw)  +  Hev(-jw)‘]/7[Hev(j«)]  (4.11c) 

Proof:  Oeflne  jJ.):  CmxB  ♦  R  by 

jj_  ( A)  *  y  M  A+A  ) 

where  *  denotes  conjugate  transpose.  Then,  using  definition  (2.8)  with  (4.10) 
-  (4.11)  we  obtain 

)ilHey^«^]  *  Ul^y^  *  fley{j“)Hey(i“l] 

>  u£Hev(j«)]  -  7(flev(ju,)RHev(jtt)]  >  0  . 

Hnece,  Heye  SPRB  . 

Comments 


(1)  In  order  to  apply  Lemma  4.1  it  is  necessary  to  have  a  detailed 

description  of  how  the  plant  uncertainty  a  propagates  onto  the 

interconnection  uncertainty  Hey  .  This  type  of  uncertainly  propagation  was 

explored  In  depth  by  Safonov  [25]  and  more  sophisticated  expressions  then 

(4.4b)  are  available  to  describe  the  uncertain  operator  Tf  .  Section  5 

ev 

contains  more  detail  on  this  issue. 


(2)  In  the  scalar  case  (4.11c)  becomes 


tc(u)  <  Re[Hev(ju.)]/|Hev(ju)| 
*  COS  *  (Hev(Ju)] 


(4.12) 


Since  c  SPRm  by  assumption,  k(u>)  Is  always  positive  for  u  e  R  ;  but 

because  of  the  cosine  function,  k(w)  <  1  .  In  Section  6  we  show  that  this 
limitation  on  the  effect  of  model  error  is  easily  violated  by  even  the  most 
benign  type  of  unmodeled  (dynamics  in  the  plant.  Methods  which  overcome  this 


/  / 


limitation  are  discussed  In  Section  7.  The  requirement  that  Mu)  <  1  also 


holds  for  any  multivariable  Hev  e  SPR 
decomposition. 


G  W 
t  ev 


W  G 
ev  r 


To  see  this  let  IT  have  the  polar 


(4.13) 


where  .  Gr  are  Hermltlan  and  Uev  Is  unitary. 

a(H  )  «  otG  )  «  o(G  )  ,  It  follows  that 
ev  t  r 


Since 


k(u»)  <  a[W  (J«)]  <  1 


(4.14) 


In  the  case  of  scalar  systems,  the  condition  k(w)  <  1  can  be  interpreted  In 

terms  of  a  limitation  on  relative  degree  of  Hev(s)  .  A  necessary  condition 

for  H  e  SPR  Is  that  the  relative  degree  of  H  (s)  does  not  exceed  one 
ev  ev 

l.e.,  phase  limited  to  ±90*.  Rohrs,  et  al .  [12]  show  that  this  necessitates 

precise  knowledge  of  plant  order,  and  hence.  Is  not  a  feasible  requirement  In 

the  presence  of  an  unstructured  uncertainty  (2.12),  where  the  order  Is 

unknown.  In  the  multivariable  case  It  Is  awkward  to  talk  about  relative 

degree  or  phase,  however,  (4.14)  expresses  the  same  limitation. 

(3)  In  several  Instances,  e.g.,  [9], [26], [27],  It  has  been  reported  that 
the  SPR  condition  has  been  eliminated.  In  each  case,  however,  it  can  be 


verified  that  the  operator  H. 


positive  constant  ,  which  is  SPR.  But, 


these  studies  do  not  account  for  unmodeled  dynamics,  thus.  In  the  notation  of 
(4.10),  only  IT  ■  positive  constant  .  Lemma  4.1  then  provides  the  means  to 
evaluate  the  effect  of  unmodeled  dynamic. 


.V'V.V 


»» 


W.V 


■' 


5.  APPLICATION  TO  MODEL  REFERENCE  ADAPTIVE  CONTROL 


Consider  the  model  reference  adaptive  control  (MRAC)  system,  shown  In 
Figure  5.1,  consisting  of:  an  uncertain  scalar  plant  P  c  RQ(s)  ;  a 

reference  model  e  Sg  ;  and  filters  with  F  e  •  The  plant  Is 

affected  by  a  disturbance  d  and  a  reference  command  r.  The  system  eauatlons 
are: 


y  -yr 

(5.1a) 

H  r 
r 

(5.1b) 

d  +  Pu 

(5.1c) 

-  e'z  a  -(8{*i  +  92z2 ^ 

(5. Id) 

F  u,  z2  *  F(y-r) 

( 5 . le) 

Assume  that  the  adaptive  law  Is  given  by  (3.15),  thus, 

ft 

a  *  r  z  e  (5. If) 


Let  the  plant  uncertainty  be  described  by(2.12),  l.e.. 


(5.1g) 


where  P*  e  R  (s)  Is  a  tuned  parametric  model  for  P.  Let  the  filter  dynamics 
*  o 

be  given  by 


F{  s) 


_  \  I 

rrsr ) 


( 5 . 1  h ) 


where  L($)  Is  a  stable  monlc  polynomial  of  degree  t  .  Thus, 

( t) ,  e2(t)  e  Ra  and  so  i(t)  e  RZl  .  Using  the  definition  of  tuned  gain 
(3.2)  we  get. 


HJi) 


v  :*  e'z  from  (3.6) 


u  »  -  e'z  «  -(e**e)'z 

•  +  ^  ~  v  • 
**1  ^*2 

*  -  -j—  u  ♦  -£—  (r-y)  -  v 


Finally, 


u 


1 

m;x/u (r_y)  "  r+Tqrr  v  c*(r_y) 


mz7U 


(5.2) 


where  A*  an(j  ^  are  p0-|ynom-( al s ,  each  of  degree  t-1  ,  whose  coefficients  are 
the  elements  of  the  tuned  gains  e*j  and  0*2  .  respectively;  and  C*  denotes 
the  tuned  controller.  The  tuned  system  (  0«O  )  Is  shown  In  Figure  5.2. 


In  terms  of  the  uncertain  plant  P,  the  adaptive  error  system  (Fig.  3.4) 
correspondl ng  to  this  MAC  system,  has  tuned  signals: 

e*  -  (1  +  PCJ’ld  ♦  [(l+PC^)-1PC*-Hr]r  (5.3a) 

(5.3b) 

and  Interconnections: 


F(l+PC*)’1C*(r-d) 

F(l+PC*)‘1(d-r) 


Hev  -  (l+PC*rlP(l+An/irl 


(5.3c) 


H  » 
zv 


F(l+PC*)“l(l+A*l/L)'1 

Fd+PCj^Pd+A^/l)*1 


(5.3d) 


The  error  system  can  also  be  described  so  as  to  highlight  the  model  error 
a  .  The  following  definitions  are  convenient: 

T.  <t*l>,C,rlP.C.  1  -  S.  (5.4.) 
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(5.4b) 


K*  Hey|  -  ( l+P*C*)"lP*( 1+A^j /L J”1 
Thus,  the  error  system  (5.3)  can  be  a) so  be  expressed  as: 

e*  -  S*(l+aT*)_1d  +  (T*(l+a)(l+aTw)"1-Hr)r  (5.5a) 


F  S*C*( 1+aT*)-1{ r-d) 
K  * 

F  S*(l+aT*)_1{d-r) 
Hev  *  K^l+aHl+aT*)"1 


35.5b) 

(5. 5c) 


F  M^d+aTJ-1 
F  K*(l+a)(l+aT*) 


(5.5d) 


The  result  that  follows  In  Lemma  5.1  gives  conditions  under  which 

Hev  e  SPRo  and  Hzv  E  SotXl  *  desP*te  mo<le^  error;  thus  conditions  (A1)-(A3) 
of  Theorems  1A  and  2B  are  satisfied.  Additional  requirements  are  necessary  to 
establish  the  class  of  tuned  signals  e*  and  z*  as  given  by  (5.5a)  and  (5.5b), 
respectively.  These  requirements  are  discussed  following  Lemma  5.1. 


Lemma  5.1:  For  the  adaptive  system  (5.3)  or  (5.5)  Hgv  e  SPRq 

H  e  S2*xl  If  the  following  conditions  are  all  satisfied: 
zv  o 


(D 


P*(s) 


g(sn’1+  Sjs""2*  ...  ♦  8n-1) 

n  fPI 

S  +ai  S  +  •  •  •  +  a* 

1  T1 


gN*(s) 


and 


(5.6a) 


(II)  Nt($)  Is  a  stable  monlc  polynomial  (5.6b) 

(III)  g  >  0  (5.6c) 

,  t  g  K,(s) 

(1v)  K*($)  ■  YTsj  c  SPRq  where  k1  ( s)  and  K?(s)  are  monlc  stable 


polynomials 


(5.6d) 

(5.6e) 


(v)  i  »  deg  L(s)  >  n  +  deg  Kj(s)  -  1 

(vl)  Ac  8^(6 )  Is  such  that 

6  (<d }  <  4(«)  :*  n(w)[n(w) |T*( ju) I  +  |S*(ju)|] 
n(u)  :*  cos  *  [ K ^ ( j«a ) ] 


-1 


V  M  c  R, 


(5.6f) 


Proof:  See  Appendix  B. 

Olscusslon 

(1)  Condition  (i)-(v)  of  Lemma  5.1  are  restatements  of  known  results, 
but  normally  they  apply  to  the  actual  plant  P,  e.g.  [7],  In  Lemma  5.1, 
however,  these  conditions  apply  to  the  parametric  model  P»  —  not  to  the 
actual  plant.  As  such,  they  are  easier  to  satisfy,  since  the  parametric  model 
Is  somewhat  arbitrary.  This  flexibility  Is  penalized  by  an  Increase  In  model 
error.  For  example.  If  the  actual  plant  has  a  relative  degree  of  2,  then 
choosing  a  parametric  model  of  relative  degree  1  —  as  required  by  condition 
(1)  —  Incrases  the  high  frequency  model  error. 

(2)  Condition  (vl)  Imposes  an  upper  bound  ?  on  the  model  error 
associated  with  the  chosen  parametric  model.  This  condition  simultaneously 
Insures  that  Hey  e  SPRQ  despite  model  error,  and  that  the  tuned  system  Is 
stable  (see  proof  in  Appendix  B). 

(3)  It  Is  easily  verified  that  3lu)  <  1  ,  as  was  discussed  following 
Lemma  4.1.  In  fact,  even  the  "optimally  tight"  bound  (see  [25]  for  details  on 
this  calculation)  given  by, 


Trrnr 


[-I1-TI  ♦  (ll+Tl2  ♦  4n  Re(KT/|K|)1/Z] 


(5.7) 


Is  also  restricted  to  be  less  than  i.  This  limitation  severely  restricts  the 
type  of  admlssable  model  error.  This  Issue  Is  pursued  in  Section  6. 


’•  • 

4Jt 


f  * 


l+i 


P. 


•'J 


I 


K 
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(4)  To  guarantee  global  stability  using  the  adaptive  law  (5. If), 
property  (1)  of  Theorem  l  requires  that  e*  ♦  0  and  z*.  e  L2*  for  all  r 
and  d.  For  example,  let  r  and  d  be  any  bounded  signals  such  that 

r  ♦  constant  and  d  ♦  constant  as  t  ♦  -  .  Property  (1)  of  Theorem  1  Is 
satisfied  If: 

6(0)  *  0  (5.8a) 

TJO)  »  H  (0)  *  1  (5.8b) 

r 

Zero  model  error  at  0C  (5.8a)  Is  certainly  to  be  expected  from  even  the  most 
crude  tuned  parametric  model . 

(5)  Let  r  be  bounded  such  that  r  ♦  constant  as  t  ♦  •  ,  but  let  d  be 
just  bounded,  l.e.,  d  e  La  .  In  this  case  It  Is  not  possible  to  guarantee 

e .  ♦  0  ,  but  we  can  guarantee  that  e*eL  .  To  obtain  global  stability  In 
this  case,  requires  the  Introduction  of  the  retardation  term  (3.16)  Into  the 
adaptive  law  (5. If),  see  part  (11)  of  Theorems  1A  or  18. 

(6)  It  Is  possible  to  obtain  versions  of  Lemma  5.1  for  adaptive  systems 

of  different  forms,  e.g..  Indirect  adaptive  [5].  Also,  the  use  of 

"multipliers",  e.g.  [4],  can  be  accounted  for  as  well.  The  multiplier 

* 

effectively  makes  use  of  the  availability  of  e  as  a  signal;  and  this  allows 
rel  deg  (P*)  *  2  rather  than  1  as  required  by  condition  (1)  of  Lenina  5.1. 


6.  LIMITATIONS  IMPOSED  BY  THE  SPR  CONDITION 


The  fact  that  the  model  error  bound  given  In  condition  (vl)  of  Lemma  5.1 
can  not  exceed  one  has  unfortunate  consequences. 

Example  1 


Consider  a  plant  with  transfer  function, 

«•>  *  p*,s)  nnrW  ,6-" 

where  P+  Is  the  parametlrc  model,  with  two  unmodeled  stable  poles  at  -a  and 
-b.  Suppose, also,  that  b  Is  much  greater  than  a,  and  that  a  Is  much  greater 
than  the  bandwidth  of  P*(s)  .  This  situation  seems  benign  —  and  most  likely 
a  certainty.  Comparing  (6.1)  with  (S.lg)  gives. 


5((d) 


.[  jl  jl/2  ,  ! 

u2*t2)d  *i>2) 


for  all  frequencies  u  >  (ab/2)*^  ,  thus,  condition  (vl)  of  Lemma  5.1  is 
violated,  and  global  stability  cannot  be  guaranteed.  The  following  example 
illustrates  this  point. 


Example  2 


Consider  the  example  MRAC  system  (Fig.  5.1)  studied  by  Rohrs  et  al .  [12], 
where: 


P(s) 


2  229 

(s+15)2  +  4 


hr(s)  •  ft 

A  A 

u  *  -  9jy  ♦  e2  r 
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o: 

to. 


► » 

kV 


v\ 

K\ 


,S 


«i  *  ye,  8^(0)  ■  .65 
e2  «  -r  e,  e2(0)  ■  1.14 

Let  r  *  constant  and  d  *  0.  Thus,  e*  ♦  0  exponentially  when  the  tuned  gains 
are  such  that  (5.8)  Is  satisfied,  l.e., 

T*(0>  *  *  Hr(0)  ■  1 

Even  though  (e*. ,  e*_)  exist  to  satisfy  this,  H  (s)  Is  not  SPR,  and  so 
l  t  ev 

global  stability  Is  not  guaranteed.  Simulation  runs  with  r  •  A  and  r  *  4.0 
are  shown  In  Figures  6.1  and  6.2,  respectively.  With  the  small  Input  (Fig. 
6.1)  we  see  a  stable  response  which  tracks  the  reference  very  well.  With  the 
large  Input  (Fig.  6.2)  the  response  Is  still  stable,  but  large  oscillations 
are  taking  place.  Larger  Inputs  will  eventually  drive  the  system  unstable, 
e.g.  [12]. 

In  this  example.  If  the  tuned  model  Is  taken  to  be  P*(s)  »  l/(s+l)  then 
It  Is  easily  verified  that  model  error  5(«)  Is  greater  than  one  at  some 
frequency. 


t 


* 


i: 
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7.  SPR  COMPENSATION 


In  this  section  we  heurlstlcally  develop  a  means  to  obtain  global  robust 
adaptive  control .  Since  the  SPR  condition  Is  violated  whenever  model  error 
exceeds  one,  a  natural  scheme  Is  to  construct  an  SPR  compensator  which 
allelvates  the  problems  by  "filtering"  the  plant  output;  thus,  avoiding  the 
trouble.  However,  direct  filtering  does  not  change  the  size  of  model  error. 
For  example,  with  the  plant  p  »  (1+a)P*  ,  let  yw  denote  the  output  of  the 
filtered  plant,  where 

y^  :*  Wy  *  Wd  +  (l+a)WP#u  (7.1) 

Thus,  model  error  Is  uneffected.  Even  filtering  Hgv  directly  by  W  offers  no 
help,  since  the  bound  (4.4c)  Is  still  less  than  one,  i.e., 

|Hey|  <  Re(W  ITev)/|W  HgJ  <  1  (7.2) 

for  any  stable  W.  What  we  seek  is  ^n  SPR  compensator  which  only  effects  the 
unmodeled  dynamics,  but  leaves  the  paramtrlc  model  Intact. 


A  compensation  scheme,  which  offers  some  promise  as  an  SPR  compensator, 

is  that  proposed  in  [22],  as  shown  In  Figure  7.1.  To  see  the  desired  result 

suppose  that  P  *  (l+a)Pm  with  &  e  B^U)  .  Then,  the  compensator  Is 

equivalent  to  a  plant  which  maps  (u,d)  Into  y  where 

c 


(7.2a) 

(7.2b) 


Thus,  whenver  s(u)  >  1  ,  select  W(s)  such  that  |w(ju)|{(u)  <  1  •  The  filter 
W  acts  like  a  "frequency  switch"  whose  function  Is  to  insure  condition  (vl)  of 
Lemma  S.l. 


There  are  two  ways  to  Implement  this  compensator  In  an  adaptive  system. 

The  first  way  Is  to  use  a  fixed  model  of  the  plant  for  Pm.  I.e.,  Pm  ■  P  . 

The  second  way  Is  to  replace  P  with  an  adaptive  observer,  I.e.,  p  •  p  . 

m  m 
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In  either  case,  to  obtain  the  benefit  of  the  SPR  compensator,  the  signal  to  be 
controlled  Is  the  compensator  output  y  ,  not  the  Diant  output  y.  Both  of 
these  compensators  will  now  be  examined. 

Fixed  SPR  Compensator 

Let  P  ■  P  ,  a  fixed  model,  and  let  the  actual  plant  be  given  by  (2.17), 
m 

P  *  (l+a)P*  with  a  e  8SU)  .  Then  the  fixed  compensator  plant  equivalent 
model  error  (7.2b)  Is: 


Pc-p* 

ac  ■■-rrc  w 


(7.3a) 


where 


iF(  jai)-P^(  jdi)  | 

SjU)  :»  |W(  jut)  1 6  (iu )  ♦  ll  -  W(jo.  11  '  I  MW  I 


(7.3b) 


This  scheme  Is  motivated  by  the  fact  that  at  low  frequencies  the  tuned 
parametric  model  P*  Is  close  to  P;  thus  6  Is  small  and  W  -  1  .  At  high 
frequencies  6  Is  large  but  (F-  -  P*)/P#  Is  small,  W  •  0  and  so  Is 

small.  Of  course  the  compensator  Is  limited  If  there  Is  large  model  error  at 
Intermediate  frequencies. 

Example  2 

Example  1  Is  modified  to  include  a  fixed  SPR  compsnator  with  W(s)  » 
l/(s+l)  and  Pts)  *  2/(s+l)  .  Simulation  results  with  the  large  step  command 
(r»4)  are  shown  In  Figure  7.2.  Comparing  these  to  Figure  6.2,  without 
compensation.  It  Is  readily  verified  that  the  Instability  tendencies  are 
eliminated.  Also,  direct  calculations  reveal  that  Hey  e  SPRQ  ,  thus  global 
stability  is  insured. 


Adaptive  SPR  Compensation 


An  adaptive  SPR  compensator,  together  with  the  adaptive  controller,  1$ 
shown  In  Figure  7.3.  The  adaptive  controller  Is  described  by. 


•T^ 


/Lc(s))  .  nc  «  deg  Lc(s) 


(7.4c) 


F^s)  -  (l/lc(s),...,s  c 


and  the  adaptive  observer  Is  described  by. 


X 

9. 


Wo  •  eo  =  y-  y 

n  -1 


(1/L  (s),  ....  s  0  /L  (s))  ,  n  *  deg  L  (s) 


(7.4d) 
(7.4d) 
(7.4f ) 


where  L  (s)  and  L  (s)  are  both  monic  and  stable.  To  generate  the  error 
o  c 

system  Interconnection  operators  associated  with  this  system,  let  9*c  and 

9  denote  the  tuned  parameters  with  respective  gain  errors,  e  and  9  ;  and 
*0  c  0 

let  v  :*  9*2  and  v  :*  9 '2  be  the  corresponding  adaptive  control  errors 
c  c  c  0  00 

(3.6).  By  analogy  with  the  procedure  used  In  Section  5  we  get. 


u 


’  C*<r-Fc'  -  ’c 


(7.5) 


where 


«  B*i  B*i 

y  *  -  j- —  d  +  (1 - £-  4)P*u  +  vQ 

0  0 


c* 


A*2/**c 

T^7TC 


(7.7) 


(7.6) 


.  .  8*2/^*o  .  gN* 
P*  1>B— /L”  U7 


(7.8) 


and  where  (A*j,  A*g)  are  polynomials  whose  coefficients  are  the  parameters  In 
6*  i  ,B*„)  are  polynomials  whose  coefficients  are  the  parameters  In 
e*Q  ;  and  N*,  p*  and  g  are  as  defined  by  (5.6a).  The  adaptive  error  model  Is 
given  below  In  terms  of  T#,  S*.  and  as  defined  In  (5.4).  In  addlton. 


define: 


R 


1  ♦  (W-l) 


The  tuned  signals  are: 

e*c  =  S*( 1+aRT#)“^R  d  +  (T*(l+&R)(l+aRT*)"l-Hr)r 
e*o  ’  O^’^I+aRT*)"^  +  0*L"1T*A(l+aRT*)“1r 


z 


*c 


FcA*2Lclp*1<*{ 1+iRT*  r 1 ( p‘Rd) 
FcS*(l+&RT*)_1(Rd-r) 


z*o 


FoA*2L^ lp*lK* ( 1+aRT* )_1 ( r'Rd ) 
FJ^I+aRT*)*1^  -  (l+a)r) 


The  Interconnections  are: 


K*{1+&R)( 1+aRT*)”1 
^*0#L“la( l+dRT#) -1 

FcP;1K*{l+dRT*)'1 

FKJl+dRKl+dRTJ*1 


-(l-W)S*(l+dRT*)'1 

l+(l-W)T#D*L“l(l+dRT*)’1 

FcA*2Lclp*l|(*(1"W>(1+4RT*) 
-F  S*(l-W)(l+dRT*)”1 


FqP^MI+aRT*)"1 
-FqK#( 1+a)( 1+aRT*)”1 


FoA*2L^P*1|<*(1'WH1+aRT*)"1 
-FqT*( 1-W)(1+a){ 1+aRT*)”1 


(7.11c) 


The  factor  (1+aRT*)"1  appears  In  all  the  terms  above.  The  transfer 

function  R  (7.9)  reduces  the  effect  of  unmodeled  dynamics;  however  not  exactly 

by  the  amount  anticipated,  vis  a  vis  (7.2).  This  is  due  to  addltonal  model 

error  introduced  by  the  adaptive  observer.  Nonetheless,  the  model  error 

attenuation  is  greater  than  with  the  fixed  SPR  compensator.  In  particular,  at 

low  frequencies  a  -  0  and  at  high  frequencies  R  «■  0  ,  since 

W  -  0  and  D*L_1  -  1  .  Without  further  testing  of  H  (7.11a)  it  is  not 
*  o  ev 

possible  to  state  that  H  e  SPR  at  Intermediate  frequencies.  Note, 

ev  o 

however,  that  the  nominal  value  of  H  is: 

ev 


K* 


K  * 
ev 


0 


-(l-W)S* 

1 


(7.12) 


which  Is  SPR0  provided  that  k*  e  SPR  and 

Re  K^(ju)  >  ^-| ( 1— W(  jw)) S*(  jw)  I ^  ,  u>  e  R  (7.13) 

Applying  (4.11)  to  (7.11a),  a  tedious  procedure,  would  give  an  upper  bound  on 

model  error  to  Insure  H  e  SPR  . 

ev  o 


44 


8.  CONCLUSIONS 


This  paper  has  presented  an  Input/ output  view  of  multivariable  adaptive 
control  for  uncertain  linear  time  Invariant  plants.  The  essence  of  the 
results  are  captured  In  Theorems  1A  and  2B  which  provide  conditions  that 
guarantee  global  stability.  Corollary  1  also  give  specific  L2  and  L^  bounds 
on  significant  signals  In  the  adaptive  control  system.  These  bounds,  for 
example,  can  be  used  to  guarantee  that  the  adaptive  system  performs  as  well  as 
a  robust  (non-adaptlve)  system  using  the  same  structure,  but  with  fixed 
gains.  By  distinguishing  between  a  tuned  system  and  a  robust  system,  we 
establish  formulae  which  can  be  used  to  restrict  the  minimum  performance 
Improvement  possible  with  the  same  control  structure. 

Although  the  stability  results  (Theorem  IA,  IB)  are  not  entirely  new  (see 
e.g.,  [7], [8]),  the  Input/output  setting  provides  the  means  to  directly 
determine  the  system  robustness  properties  with  respect  to  model  error.  The 
type  of  model  error  examined  can  arise  from  a  variety  of  causes,  such  as 
unmodeled  dynamics  and  reduced  order  modeling.  It  Is  very  difficult  to  treat 
this  typ«  of  "unstructured"  dynamic  model  error  by  using  Lyapunov  theory, 
since  the  system  order  may  not  be  known  —  In  fact,  it  may  be  infinite. 
Although  infinite  dimensional  (distributed)  systems  were  not  considered  here. 
Theorem  1  can  be  modified  to  Include  them,  e.g.,  [26]. 

The  structure  of  Theorems  1A  and  IB  require  that  a  particular  subsystem 
operator,  denoted  Hey  ,  Is  strictly  positive  real  (SPR).  This  requirement  Is 
not  unique  to  this  presentation  -  passivity  requirements.  In  one  form  or 
another,  dominate  proofs  of  global  stability  for  practically  all  adaptive 
control  systems.  Including  recursive  Identification  algorithms. 

Unfortunately,  although  Hey  e  SPR  is  robust  to  model  error  (Lemma  4.1),  the 
bound  on  the  model  error  Is  too  small  to  be  of  practical  use.  Even  the  most 
benign  neglected  dynamics  violate  the  bound. 

Although  this  paper  is  concerned  with  continuous- time  systems,  the 
theorems  carry  over  virtually  Intact  to  discrete-time  systems.  This  Is  a 
direct  consequence  of  the  portable  nature  of  the  Input/output  view.  However, 
there  Is  an  Important  Issue  unique  to  discrete- time  systems:  plant 


uncertainty  is  critical  to  where  performance  is  actually  measured,  which  is  in 
continuous- time,  not  at  the  sampled-data  points.  As  a  consequence,  it  may  be 
necessary  to  map  the  discrete  portions  of  the  adaptive  system  (most  likely  the 
controller)  Into  continuous-time,  l.e.,  the  L2"9*(ns  of  the  discrete-time 
operators  In  the  Interconnection  map,  which  are  associated  with  the  adaptive 
discrete-time  controller,  would  be  needed  rather  the  discrete-time  ig-galns  • 

Another  area  worth  pursuing  is  the  adaptive  control  of  non-linear 
plants.  The  plant  uncertainty  description  (2.11)  does  not  exclude  non-linear 
plants.  Mote  that  slowly  drifting  parameters  in  an  otherwise  perfectly  known 
LTI  plant  could  yield  the  same  uncertainty  description  as  a  non-linear  plant 
approximated  by  a  parametric  LTI  model.  All  that  is  required  is  that  there 
exists  a  (possibly)  Infinite  dimensional  LTI  system  which  matches  the 
input/output  behavior  of  the  plant  for  each  possible  input/output  pair.  Of 
course,  if  the  plant  Is  truly  non-linear,  then  the  tuned  control  Is  likely  to 
be  non-Hnear,  which  raises  some  very  interesting  issues  for  further  research. 

One  final  remark:  the  stability  results  presented  here,  as  well  as  other 

known  results,  provide  global  stability.  This  Is  achieved  by  requiring 

H  e  SPR  ,  a  condition  which  is  difficult  to  maintain  In  normal 
ev 

circumstances.  On  the  other  hand,  this  Is  a  sufficient  condlton;  violation  of 
which  does  not  necessarily  lead  to  Instability.  The  simple  example  presented 
here  In  Figure  6. 1-6. 2,  Illustrates  the  point.  Other  examples  of  this 
phenomena  abound,  e.g.,  [12].  It  would  appear  then,  that  a  more  valid 
approach  to  providing  a  system- theoretic  setting  for  adaptive  control  Is  to 
develop  local  stability  conditions,  which,  hopefully,  do  not  require  that 
Hev  e  SPR  .  Preliminary  results  on  local  stability  support  this  hope,  e.g., 
[33],  [34]. 


ACKNOWLEDGEMENT 


We  wish  to  acknowledge  many  useful  discussions  on  portions  of  this 
material  with  C.R.  Johnson,  Jr.  (Cornell  University),  8.0.0.  Anderson 
(Australian  National  University),  C.E.  Rohrs  (Notre  Dame  University),  and  C.A 
Desoer  (University  of  California,  Berkeley). 


APPENDIX  A 


PROOF  OF  THEOREMS  1  AND  2 

Prel Imlnarles 

The  main  Ingredient  In  the  proof  Is  to  show  stability  by  means  of 
passivity.  Although  there  are  many  variations  on  this  theme,  a  general  result 
Is  given  by  the  following. 

Theorem  A.l  ([21],  [35] 

Consider  the  feedback  system  of  Figure  A.l  below  with  causal  operators 
Gj  and  G2  • 


Figure  A.l  Feedback  System 


Suppose  there  exists  real  constants  e^,  <5^ ,  ,  1*1,2  ,  such  that 

<x,GjX>t  >  €{1X1*2  +  «{lG|X»^2  +  a{,  Vt  >  0,  V  x  e  LgCO.tl 

for  1*1,2.  Then  the  following  holds  V  t  >  0, 

2  2 

(cg+djiiyji^  ♦  (c1+«2J|y2,t2  <  ,yl,t2(,ul,t2  +  ** ,u2at2) 

2  2 

♦  iyzit2(iu2it2  +  2|c1l«iu1it2)  ♦  |c1I*iu1i ^  ♦  lc2l  •  i 

+  |o,|  ♦  la- | 


(A.l) 


(A. 2) 


Proofs  of  both  theorems  also  rely  on  well  known  results  for  systems 

H  e  Snxm  .  The  results  required  here  are  summarized  In  the  following, 
o 

Theorem  A-2  [see  [19],  Thm.  9,  pg.  59] 

Let  H  c  Snxm  ;  then: 

0 

(1)  If  Utlj  ,  then  y  -  Hu  C  Lj  L^  ,  *  e  .  y  Is  continuous,  and 
y(t)  ♦  0  as  t  ♦  •  . 

( il )  If  u  c  Lm  ,  then  y  ■  Hu  e  -n  .  9  c  Ln  .  ««d  Y  Is  uniformly 

mm  urn  mm 

continuous. 

(Ill)  If  u  e  L™  and  u(t)  ♦  constant  c  e  Rra  as  t  ♦  -  ,  then 

y(t)  ♦  H(0)c  exponentially  as  t  ♦  •  . 

In  order  to  simplify  notation  we  drop  the  superstrlct  on  l”  which 
Indicates  vector  size. 

We  will  establish  Theorem  1A  first.  Some  of  the  steps  will  be  repeated 
for  IB.  Also,  without  loss  of  generality,  the  matrix  r  in  the  adaptation  law 
(3.151,(3.16)  Is  set  to  Identity.  Corollary  1  is  established  as  a  by-product. 

Proof  of  Theorem  1A 

Part  (1) 

Identify  G^,  Gg  In  Figure  A. I  with  e  ♦  v  and  Hey  respectively.  Also, 

let  u  -  V  u2  •  0,  et  «  e,  y1  ■  e2  ■  v,  and  y2  ■  H#yv.  .  Using  adaptive 

law  (3.15)  we  obtain, 

<e,v>T  ■  <e,Z'e>T  ■  <Ze,o>T  ■  <e,  e>T  (A. 4) 

*  ■y  i9(T)i^  -  y  iq(0)|2  (A. 5) 

)  ■  y I0(O)|^  (A. 6) 


fH 


I  •  • 

r. 


Thus,  using  (A.l)  gives. 


1  2 

el  *  *1  *  0,  “1  *  "  7  ,e{°”  (A. 7) 

Since  G0  ■  H  e  SPR.  by  assumption,  3u,  y  >  0  such  that  1  xe  L,  , 

C  CV  y  ▼ 

<x,Heyx>T  >  iii xi ^2 ,  |Hevx,T2  *  YlX|T2  *  H*nce*  from 


€2  *  u’  52  *  °2  *  ®  (A. 8) 

Using  Lemma  A.l,  together  with  (A.4MA.8)  gives, 

l*lT2  <  *  C * *** T2  *  2«la<0)|2)l/Z]  (A.9) 

ie-e*iT2  <  Yiviyg  (A. 10) 

le (T) 1 2  <  le (0) | 2  ♦  2ieiT2  ivit2  (A. 11) 


The  bounds  shown  In  (4.8)  follow  using  the  assumption  e*  e  L2  .  Hence, 

e,v  c  L.  and  9  e  L  . 

z  • 

.  z 

Having  established  that  vcL,  ,  Theorem  A-2  — ►  z:»  z- z*  e  L0Ol  ,  z  e 

,  .  c  c.  ■ 

z  ♦  0,  and  z  Is  continuous.  Since  z*,  z*  e  L  by  assumption,  it  follows 
that  z  t  and  i  e  L  (-■►  z  Is  uniformly  continuous).  Using  v  *  Z'e  with 
z,  81  el  4veI  .  Using  e  «  e*«H  v  with  e*  e  L  and  H  e  S  (by 
assumption),  and  vc  1^4  et  .  Hence,  fi  •  Ze  e  9  *s  uniformly 

continuous^  v  »  Z'e  is  uniformly  continuous  (since  z  1s)^  v  ♦  0  since 
v  e  Lg  Is  established.  Using  v  ♦  0^  e  -  e#  ♦  0  ,  and  since  e*  ♦  0  by 

assumption,  e  ♦  0  .  Furthermore,  v  ♦  ()■*  z  ♦  0  exp.  and 

a  ■  Ze  ■  Ze  +  Z*e  ♦  0  ,  because  z  and  e  ♦  0  .  Using  v  »  Z'e  +  Z'§  with 

z.  9,  a  e  L^-4  V  e  Lw  .  Hence,  e*  «  6*  -  Hev  1  c  ,  because  e*  e  L_  by 

assumption.  Thus,  e  «  2e  +  Ze  e  l  •  This  establishes  properties  (1-a)- 
( 1-d) . 

To  show  (1-e)  consider  (3.15)  written  as: 


(A. 12) 


^  ~  »  A  -  V  "  k  •-  *  ■ »  ^  ^  i*  \.  «.  -  ^  .* 


®  ■  -Z.  Ney  K  >  ♦  » 


»:*  -(Z."„Z‘  ♦  Z  Hev  z;  ♦  Z  He>  Z')6 

Since  we  have  already  established  that  z  *  0  exp.  and  e  e  L  .  It  follows 
that  w  ♦  0  exp.  Since  z*  e  PE  by  assumption  (provided  e*  »  0)  ,  w  Is  exp. 
stable  by  (2.15).  Hence,  e -*•  0  exp.  e,  v  ♦  0  exp.-^  e-e*  ♦  0  exp.  This 
completes  the  proof  of  part  (1)  with  adaptive  law  (3.15). 

To  show  that  (l-a)-(l-a)  hold  with  adaptive  law  (3.16)  requires  showing 
that  Gj:e  ->  v  Is  passive.  Consider  the  typical  time  Interval, 


I 


Ij  *  { t  e  CtQ.tj) j  ie( t)i  <  c} 

I2  *  {t  e  Ctj.tg)!  ie(t)i  >  c  >  maxie*i} 


Hence, 


Thus, 


<e,v>  *  <e,v>  +  <e,v>T 

i  lt  i2 


(A. 13) 


(A. 14) 


<e,v>.  -  <e 

ll 


<e,v>. 

l2 

■  j  1 9 ( t2  )■ 2 
>  j  leftj)!2 


•  0>I1  *  7  "  7  ,0(to)|2 

•  *  2  * 

■  <8  ♦  (1  -  iei/c)  e,  e>T 

l2 

-  ^  i a( t^ )i 2  ♦  (1  -  i0i/c)2<e,  e>j 

-  y  ie( tj )■ 2 


2 


(A. 15) 
(A. 16) 
(A. 17) 
(A. 18) 


because  <e,  e>T  >  0  from. 


r  »  -  TTv'  i.~H  iTW  i.TE  "W  -"V  ■"%»  I. 


e( t) *  e(t)  «  0 ( t) '[e(t)-e#] 

*  ■ 0 ( t) i 2  -  8(t)'a* 

>  ie(t)i2  -  ie(t)ic 

*  ie(t)l (l8<t)i>  c)  >0,  Vt  e  Ij  •  (A. 19) 

Thus, 

<e,v>j  >  j  i 8 ( t2 )■ 2  -  j  i 8 ( tQ ) i 2  (A. 20) 

Repeating  the  above  procedure  recurslvely.we  eventually  conclude  that 
1  2 

<e,v>T  >  -  y  i 8 ( 0 ) i  as  before  {A. 6),  and  hence,  Gfe  h»*v  Is  passive.  The 
results  In  (1)  now  repeat  for  adaptive  law  (3.16).  This  completes  the  proof 
of  part(1). 

Proof  of  Theorem  lA,Part  (11) 

Theorem  1A,  Part  (11)  Is  essentially  an  L^-stabllity  result.  The  method 
of  proof  requires  the  notion  of  "exponential  weighting"  which  Is  a  means  to 
obtain  L^-stablllty  of  a  system  from  the  l^-stablllty  of  an  exponentially 
weighted  version  of  the  system  (see  e.g.,  [19],  Chapter  9).  We  require  the 
following: 

Definition;  Given  a  real  number  a  define  the  exponential  weighting  operator 
by 

xa(t)  :«  eotx( t)  (A. 21) 

Consider  the  system  y  ■  Gu.  An  exp.  weighted  version  of  this  system  Is 
denoted  by  y°  :■  G®  u°  .  Mote  that  If  G  is  a  convolution  operator  with 
transfer  function  G(s)  then  Ga  Is  also  a  convolution  operator  with  transfer 
function  G(s>a)  .  Thus,  the  correspondl ng  exponential  weighted  error  system 
corresponding  Is  described  by 
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(A. 22) 


*  *  * 


2«  .  2«  .  H«v  v® 


v®  *  z'e® 

a®  *  08®  ♦  Ze®  -  p(a)e® 


where  a  >  0  such  that 


H®  e  SPR™  and  H®  e  S** 
ev  +  zv  o 


(A. 23) 


Using  Theorem  A-l,  Identify  with  e®  ♦  v®  and  with  H^.  Mote  that  It  Is 
always  possible  to  find  some  a  >  0  such  that  (A. 23)  holds.  We  now  examine 
the  passivity  of  G,:  e®  ♦  v®  .  Thus, 


<e  ,  v®>,  «  <e  ,  Z‘a®>T  *  <  Ze®.  a®>T 

a  *a  a  , *.*a 
*<9,0  -  o0  +  p(0)0  >j 

*  y  I0(T)|  -  y  I0(O)|  ♦  <p(0)0a,  0a>j-ol0°ly2 

>  y  e^®^l0(T)l^  -  y  1 9 ( 0 ) | ^  -  alO®lTp 


(A. 24) 


The  last  line  follows  from  (A. 19),  hence,  (A. 24)  holds  with  or  without  the 

retardation  term  In  the  adaptive  law.  At  this  point  there  are  two 

possibilities:  either  8  t  or  |8(t)|  ♦  •  as  t  ♦  •  •  If  9  el.  then  ; 

constant  c  <  •  such  that  191  <  c  .  Then, 

0  «o 


<e®,  v“>T  >  j  e2oT(i8(T)i£  cj)  -  £  10(0)1 

1  2aT  2  1  ..mi.2 

>  -  j  c  co  "  7  19(0)1 


2  1 


(A. 25) 


If  | a ( t ) |  ♦  •  as  t  ♦  •  then  It  Is  always  possible  to  select  an  arbitrarily 
large  T  such  that  ie(T)i*  I8IT-  .  Hence,  for  this  T,  (A. 24)  becomes, 

<e®,  v®>.  >  i  e2aT(i9(T)i2  -  *0iyj  -  i  «9(0)« 2 


I 


Thus, for  some  arbitrarily  large  T,  (A. 25)  and  (A. 26)  have  the  general  form 
l.e.. 


<e°,  v“>T  >  -c1  e2oT  -  c2  -c(aT) 


where  Cj,  c2  are  non-negative  constants.  Hence, 


el  *  Sl’  *  °»  ai  *  -c(aT) 


Since  G2  *  H^y  c  SPR+,  B  constants  u,  y  >  0  such  that 


<x.  Hjy  x>T  >  uix,‘2 

,Hev  X,T2  <  7  ,X,T2 


Then, 


*  u .  62  *  «2  *  0 


Using  (A. 2),  we  get 


,ya,T2  <-7i,eS,T2  +  (,e2,T2  +  2u  C(«T))1/2] 


Since  e*  e  L  by  assumption, 

i e*i j2  «  e°T(2«)"1/2ie*i- 

Thus, 


-  £lU»Tl^_r..  .  w  —  .2  .  ,  ,-Zc.T tM1/2i 


IV-,T2  <  -  .-£■ - [le*^  ♦  (ie*r  +  4a  e'^ydaT))1'4] 


Since  H^v  c  SgXm  ,  we  obtain 


l*(T)l  »  I  £  H2v(T-T)v(T)dr| 

-  le"*1  ^¥(T-T)va(T)dT| 


(A. 27] 


(A. 28] 


(A. 29) 


(A. 30) 


(A. 31) 


(A. 32) 


(A. 33) 


(A. 34) 


{A. 35) 


(A. 36) 


<  ®*'8,Ti»£v(-)i1  •  »va« 


T2 


where  H^y(t)  Is  the  Impulse  response  matrix  associated  with  H^y 
Substituting  (A. 33)  and  (A. 27)  Into  {A. 36)  and  noting  that 
e"^  c(aT)  <  Cj  +  c2  ,  we  obtain. 


|z(T) I  <’(2a>~-— ■  iH^v(.)it  .  [i e*iw  ♦  (ie*«2  ♦  4«u(c1+c2))1/2] 


(A. 37) 


Since  the  right  hand  side  Is  Independent  of  T,  and  since  T  can  be  selected  to 
be  arbitrarily  large.  It  follows  that  i  t  L<  .  Assuming  there  Is  no 
retardation  or  persistent  excitation,  this  completes  the  proof  of  (11-a)  to 
(11-d). 

Assume  now  that  z  c  PC  ,  which  Is  a  noncontradictory  assumption  since  we 
have  already  shown  that  z  e  L  .  Hence, 


§  -  -Z  Hey  V  e  ♦  Ze* 


(A. 38) 


Since  z  e  PE,  H  e  SPIT  and  z,  e*  e  l  ,  It  follows  from  (2.15)  that 
(Ze*.  9(0))  (-►0  Is  exp.  stable,  thus,  a,  5  e  L#  .  The  remaining  results  In 
(11-e)  follow  Immediately. 

Suppose  now  that  the  adaptive  law  Is  given  by  (3.16).  Then,  we  can 
wrl te. 


e  *  Z  e  -  p(e)e  ■  Z[e*-HevZ'(9-e*)]  -  p(e)e 


w  -  z  Heyz'  e  -  p(e)e 


(A. 39) 


where  w  :*  Z  e  ♦  Z  H  Z'9*  e  l  ,  because  z,  e*  c  L  •  Consider  the 

ev  w  •  •* 


candidate  Lyapunov  function  V:tH^i9(t)i  .  Hence, 


V  -  2  w'e  -  o'  Z  HeyZ'9  -  p( 0 )V 


(A. 40) 


Suppose  ie(t)i  ♦  •  as  t  ♦  •  .  Then  there  exists  a  time  T  >  0  such  that 
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ie(T)i  ■  I9«T-  *  VT  >  c  .  Hence, 

VT  <  aiwi^  vJ/2+  iz»2  y.(h6V)vt  -  (i  -  yJ/2/c)2vt  (a. 4i) 

Clearly,  there  exists  a  finite  constant  c^  such  that  when  V  >  c^,  <  0  . 

Therefore,  a  can  not  grow  beyond  all  bounds,  and  hence,  e  E  ^  .  So  then  Is 
e  and  3  ,  and  again  the  result  of  (11 -e)  follow.  This  completes  the  proof  of 
Theorem  1A.  Note  that  In  this  case  we  do  not  obtain  specific  bounds  on  e, 
because  the  proof  proceeds  by  contradltlon. 

Proof  of  Theorem  IB 

Part  (1) 

Since  H  c  SPR  ,  there  exists  q  >  0  such  that  G  :*  {1  +  qs)H  e  SP 
ev  .o  ev  ev 

and  furthermore,  G“*  e  S  .  As  a  result  we  can  write  (3.14a)  as, 

e  *  -Hev  y.  y  •  »  -  o;J(e.  *  q  *.)  <»•«) 

Referring  to  Lenina  A-l,  let  Gj  :  vi->e,  G2  *  Hev  ’  ui  *  0  *  and 
u2  *  -G”*(e*  +  qe*)  .  Using  (A. 2)  together  with  (A. 42)  and  the  passivity 
properties  of  H  gives, 

■  ®»T2  *  ^j{iu2It2  *  (,U2*T2  +  |s(0)  |2}^2]  (A. 43) 

le(T)|  <  !e(0)|  +  ZieiT2  •  1  u2* T2  (a. 44) 

where  u  Is  defined  In  (4.9a).  Using  (4.9b)  gives, 
i u2i Tg  <  (l/k)ie*  +  qe*iT2  *  Th1s  t09ether  w1th  (A. 43),  (A. 44)  and  the 
assumption  e*,  e*  c  l2  gives  the  bounds  shown  In  (4.9).  Hence, 
e  e  L,,  e  e  L  .  However,  we  can  not  conclude  that  v  e  L.  as  In  Theorem  1A, 
part  (1).  From  (A. 42),  we  can  conclude  that  (l  +  qs)  v  e  L2  •  Since 
G  :*  (1  +  qs)H  e  S  ,  It  follows  from  Lemma  A-2  that 

-ZV  ZV  Ot 

z  :■  z-z*  e  L2  n  Lw  #  z  c  L2  and  z  ♦  0  •  Repeated  use  of  Lemma  A-2  and  the 
error  equations  (3.14)  gives  the  results  (1-a)  -  (1-d).  (1-e)  follows  from 

the  arguments  In  the  proof  of  Theorem  1A,  part  (1). 
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APPENDIX  B 
PROOF  OF  LEMMA  5.1 

The  proof  utilizes  the  following  known  results: 

Definition:  Let  J  denote  a  subset  of  S,  consisting  of  functions  In  S  whose 
Inverse  Is  also  In  S. 

Fact  [29]:  If  G  is  any  scalar  transfer  function  In  R(s),  then  G  has  a  coprime 
factorization  In  S,  i.e.,  there  exists  N,  0,  A,  and  B  In  S  such  that 
G  *  N/D  and  AN  +  BO  *  1. 

Lemma  B-l:  Consider  the  tuned  adaptive  system  of  Figure  5.2.  Let 
P*  e  Rq( s)  and  c*  e  R0( s )  have  coprime  factorizations  in  S  given  by 
P*  *  Np/Dp  and  C*  *  Nc/Dc  »  respectively.  Then,  the  elements  of  the 
transfer  matrix  from  (r.d)  into  (e*.  z#,y,  u)  all  belong  to  S,  If: 

(i)  Q  :*  Dp0c  ♦  NpNc  e  J  ,  (from  [29])  (B.l) 

and 

(ID  6(<it)|T*(j<it)|  <  1,  V  u  £  R  ,  (from  [16]) 

where 

T*  :«  NpNc/Q  :-  P*C*( l+P*C*)-1  (B.2) 

Using  the  definition  of  Q  we  can  write  H  and  H  from  (5.5)  as, 

ev  zv 

Hev  "  NpQ^d+aHl+aT*)"1  (8.3) 

FOpQ-^l+aT*)-1 

FNO^d+AJd+AT*)*1 

p 


(B.4) 


From  the  definition  of  K*  (5.4b),  we  also  obtain 


r 

Q  »  MpK;1  (B.3)  r 

<e 


Proof  of  Lemma  5.1 


Me  first  show  that  (1),  (11),  and  ( 1 v)  *>  Q  e  0  .  Let  P  *  N  /D  be  a 

*  p  p 

coprime  factorization  of  P*  such  that  rel  deg  0p($)  *  0*  Since  (1)  *>  rel  deg 
P*(s)  *  1,  It  follows  that  rel  deg  Np(s)  ■  1.  Moreover,  (iv)  *> 
rel  deg  K*(s)  «  1,  and  that  <i(s)  and  are  stable.  This,  together  with 

(11)  and  (B.3)  establishes  that  Q  e  J  . 

H  e  S  follows  Immediately  by  Inspection  of  (B.2),  since:  F  e  S  by 

ZV  0  0 

assumption;  0^,  e  S  ;  Q  e  J;  a  e  S  by  assumption  (vl);  and  finally  (vl) 

*>  (11)  of  Lemma  B-l  *>  (1+aT*)”1  e  S  . 

Conditions  (Iv)  and  (vl)  »>  H  c  SPR  .  This  follows  from  Lemma  4.1 

ev  o  . 

by  letting  IT  *  and  letting  l  +  H  ■  (l+aMl+aT^)”1  .  Thus,  (4.4a)  Is 

satisfied  since  K,  t  SPR  from  (Iv).  Also,  from  (4.4b), 

*  0 

k(w)  *  |flev(j<u)|  *  la(  ju)S*(  jo>)[  l-a(  ju)T«(  jail]  M  (B.4) 


5(w)  lS*(  Jti» )  1 

<  rTCTTTTTOT  <  ^m)  "nU' 


(B.5) 


The  last  Inequality  comes  from  conditions  (vl)  and  the  definition  of 
T(u)  from  (4.4b). 


The  final  step  In  the  proof  of  Lemma  5.1  Is  to  show  that  there  are  a 
sufficient  number  of  parameters  In  e*  to  Insure  a  solution  exists.  This  Is 
guaranteed  by  satisfaction  of  condition  (v).  To  see  this  combine  (B.3)  with 
the  definition  of  Q  from  (B.i)  to  get 


Mn  +0n°r 

c  p  pc 


W*1 

P 


r  1 
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.%  /.  '•  *v 


(B.6) 


fi 


From  (5.2),  let  Mc  ■  A^A  and  Dc  »1  +  A#^  be  a  coprime  factorization  of 
C*,  and  let  Np  ■  g  M*/l  and  0p  ■  1  ♦  D*/l  be  a  coprime  factorization  of  P*» 
where  P*  is  as  defined  In  (1).  with  K*  given  by  (1v),  (B.6)  becomes  the 
polynomial  equation, 

A*!  +  **2KiH*  “  UK^-KjD^)  (B.7) 

Since  deg(K2M#)  «  degfl^Dj  and  KJt  K2,  N*.  and  D*  are  all  monlc,  it  follows 
that  deg[L(K2M#-K *  deg(L)  +  deg(K^)  +  deg(D#)  -  1  .  Then,  using  known 
results  on  polynomial  equations,  e.g.  [30],  It  can  be  shown  that  (v)  Implies 
that  (B.7)  has  a  solution  (A*j,  A*2)  . 
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An  efficient  algorithm  for  output-error  model  reduction^ 

BOAZ  PORAT J  §  and  BENJAMIN  FRIEDLANDER* 

A  new  algorithm  ia  presented  (or  reduced -order  modelling  of  linear  discrete-time 
systems,  using  an  output-error  criterion.  A  closed  form  expression  is  developed  for 
the  gradient  of  the  cost  function  with  respect  to  the  model  parameters.  A  computa¬ 
tionally  efficient  algorithm  for  computing  this  gradient  is  derived.  A  Fletcher-Powell 
optimization  procedure  utilizes  the  gradient  vector  to  compute  the  reduced-order 
model  parameters.  A  special  initialization  procedure  ia  proposed,  and  the  stability  of 
the  reduced-order  system  is  monitored.  The  performance  of  the  algorithm  is 
illustrated  by  some  numerical  examples. 

1.  Introduction 

The  problem  of  mathematical  modelling  of  physical  phenomena  arises  in 
many  scientific  disciplines.  An  important  aspect  of  modelling  is  the  conversion 
of  complex  models  into  simpler  ones.  It  is  usually  desirable  to  use  models  that 
are  as  simple  as  possible  yet  still  capable  of  capturing  the  salient  features  of  the 
underlying  phenomena.  Model  simplification  leads  to  savings  in  computa¬ 
tional  requirements  and  hardware  costs  and  facilitates  the  analysis  and  under¬ 
standing  of  complex  problems.  In  this  paper  we  consider  a  technique  for  the 
reduced-order  modelling  of  linear  discrete-time  systems. 

The  problem  of  reduced-order  modelling  (sometimes  called  rational  approxi¬ 
mation  on  the  unit  circle)  can  be  defined  as  follows :  let  g°{z)  be  a  rational 
.Vth-order  transfer  function 

J  v  '  fflO(z)  l-HV»z-‘+  ...  +«„•*-'• 

where  the  polynomial  a(z)  is  assumed  to  be  stable,  i.e.  to  have  all  its  roots 
strictly  inside  the  unit  circle.  Let 

kz),*£>-  (2) 

'  a(z)  1+0,2-'+  ...  +O.Z-"  '  ' 

be  an  nth-order  approximation  to  g°(z)  (where  n  <*V),  in  the  sense  that  g(z)  is 
‘  close  ’  to  g°(z)  under  some  criterion.  The  criterion  used  in  this  paper  is  the 
L2  norm  of  the  difference  on  the  unit  circle,  i.e. 

1  • 

J  (e*P  (j«))-0  (exp  (jtu))|*d<u 

=  J_  l  t>°  (exp  (jut))  b  (exp  (jw))  * 

2rr  a9  (exp  (jut))  a  (exp  ( jut )) 
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Using  Parseval’s  theorem,  the  criterion  can  be  alternatively  specified  in  the 
time  domain.  Let 

S  9°  :  £7(z)  =“  S  (4) 

1-1  1-1 

Then  * 

F~  I  (0,#-< 7,)*  (5) 

«-l 

Minimization  of  F  as  defined  by  (3)  or  (5)  over  all  possible  parameter  values 
{6„  at:  n},  will  determine  the  optimal  reduced-order  model  g(z).  We 

will  refer  to  this  procedure  as  the  output-error  method. 


Figure  I .  A  model  tor  the  output  error  method. 


The  name  ‘  output-error  ’  comes  from  the  system  identification  literature 
(Landau  1979).  Consider  the  problem  depicted  in  Fig.  1  :  two  systems  g°(z) 
(the  real  system)  and  g(z)  (the  model  to  be  estimated)  with  a  common  input 
process  vt,  a  unit-variance  white  noise  process.  It  is  desired  to  estimate  the 
parameters  of  <7(2)  so  that  the  mean-square  error  between  the  outputs  of  the  two 
systems  E{e ,2}  will  be  minimized.  It  is  a  straightforward  matter  to  check  that 
the  mean-square  error  criterion  is  identical  to  V  as  defined  earlier. 

The  output-error  criterion  seems  to  be  a  good  candidate  in  many  applica¬ 
tions.  It  uses  a  physically  meaningful  error  criterion  and  leads  to  satisfactorv 
performance  in  the  context  of  estimation  and  control  problems.  The  main 
difficulties  with  this  method  are  related  to  the  computation  of  the  reduced-order 
model.  First,  the  error  function  is  a  non-quadratic  function  of  the  model 
parameters  (ffl(z)).  Therefore  the  minimization  of  this  function  involves  a  non¬ 
linear  optimization  procedure.  Such  procedures  are  often  complicated  and 
computationally  expensive,  especially  for  high-order  systems.  Second,  the 
error  function  V  will  generally  have  multiple  local  minima,  making  it  difficult 
to  reach  the  global  minimum. 

A  number  of  model-reduction  algorithms  based  on  the  output  error  have 
been  proposed,  mainly  in  the  context  of  filter  design.  Sanathanan  and  Koerner 
(1963)  have  proposed  an  iterative  procedure  in  which  a  conditional  output  error  is 
minimized  at  each  stage,  where  the  conditioning  is  on  the  denominator  poly¬ 
nomial  computed  at  the  previous  stage.  Steiglitz  and  McBride  (1965)  used  a 
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similar  idea  in  a  system  identification  context.  Deczky  (1972)  proposed  a 
technique  for  minimizing  the  p-norm  of  magnitude  error 

f  (I 9  (exp  (>»))l  -  \9°  (exp  (jut)) |)*»  dw 

and  the  p-norm  of  the  phase  error 

J  |phase  (g  (exp  (jut)))  -  phase  (g°  (exp  (jo>))))p  dot 

The  reduced  transfer  function  g(z)  is  modelled  as  a  cascade  connection  of 
second-order  filters.  This  procedure  is  extensively  used  for  filter  design. 
Aplevitch  (1973)  gave  a  gradient  algorithm  based  on  a  state-space  formulation. 
Recently,  Yahagi  (1981)  proposed  a  gradient  algorithm  for  minimizing  the 
output  error  with  respect  to  a  model  specified  by  a  finite  number  of  impulse 
response  terms. 

A  number  of  alternative  procedures  have  been  proposed  in  the  literature, 
apparently  stimulated  by  the  difficulties  in  computing  the  output-error 
reduced-order  model  parameters.  Perhaps  the  most  popular  of  these  is  the 
so-called  equation-error  method  which  uses  an  error  function  of  the  form 

1  * 

V “jjj  f  1°  (*XP  0<"))  9 0  (exp  ( jut))-b  (exp  ( jut))\*dw 

]  |«  (exp  (ju>))|*|g°  (exp  (jot))-g  (exp  (jo>))|*d<«>  (6) 

Note  that  this  cost  function  involves  a  filtered  version  of  the  output  error.  The 
equation-error  method  has  the  advantage  of  being  quadratic  in  both  the  a(z) 
and  the  6(2)  coefficients  ;  hence  the  minimization  procedure  is  fairly  straight¬ 
forward.  On  the  other  hand,  the  error  function  tends  to  put  a  small  weight  on 
frequencies  where  the  magnitude  of  the  response  is  large,  yielding  poor 
approximations  for  systems  with  poles  near  the  unit  circle.  Even  more 
problematic  is  the  fact  that  a(z)  resulting  from  minimizing  V  is  not  guaranteed 
to  be  stable. 

Many  model-reduction  methods  that  are  not  based  on  the  output-error 
technique  have  been  proposed  in  the  literature.  YVe  mention  in  particular  the 
relatively  recent  development  of  the  balanced  realization  method  (Moore  1978) 
and  the  optimal  Hankel-norm  method  (Kung  1980).  Other  well-established 
techniques  include :  dominant  mode  approximation,  aggregation,  singular 
perturbation,  Routh  approximation  and  Pad6  approximation.  Here  we 
consider  only  the  output-error  method,  which  appears  to  work  well  in  various 
control  and  signal  processing  applications. 

In  this  paper  we  present  a  new  algorithm  for  the  direct  minimization  of  the 
output  error  function  (3),  (5).  Through  a  detailed  analysis  of  the  error  function 
V  we  were  able  to  develop  a  closed-form  expression  for  the  gradient  vector  (i.e. 
the  derivatives  of  V(g)  with  respect  to  the  parameters  o„  b().  This  gradient 
vector  is  then  used  in  a  Fletcher-Fowell  minimization  algorithm  to  compute  the 
parameters  of  the  reduced-order  model.  Using  some  facts  from  the  theory  of 
discrete  Lyapunov  equations  and  Toeplitz  matrices  we  were  able  to  develop  an 
efficient  algorithm  for  computing  the  gradient  vector,  requiring  of  the  order  of 
•V*  multiplications  and  additions.  This  seems  to  be  by  far  more  efficient  than 
any  other  existing  schemes  for  computing  the  gradient.  A  special  initialization 
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procedure  based  on  some  properties  of  orthogonal  polynomials  is  proposed. 
This  procedure  seems  to  provide  a  good  starting  point  for  the  subsequent 
minimization  algorithm,  leading  with  high  probability  to  the  global  minimum. 
A  unique  feature  of  the  complete  algorithm  is  that  it  guarantees  stability  of  the 
reduced-order  model  (i.e.  of  a(z))  at  each  step. 

We  believe  that  the  technique  proposed  in  this  paper  provides  for  the  first 
time  a  satisfactory  solution  to  the  output-error  reduced-order  modelling 
problem.  The  algorithm  has  a  number  of  properties  that  distinguish  it  from 
previous  attempts  in  this  direction  :  (i)  exact  closed-form  computation  of  the 
gradient ;  (ii)  computational  efficiency  ;  (iii)  improved  initialization ;  and 
(iv)  guaranteed  stability  of  the  model.  These  features  make  it  a  viable  and 
practically  implementable  technique.  Our  limited  computational  experience 
with  the  algorithm  has  been  very  favourable. 

The  outline  of  the  paper  is  as  follows.  In  §  2  we  derive  the  closed-form 
expression  for  the  gradient.  In  §  3  we  discuss  the  implementation  of  the 
method,  in  particular  the  efficient  computation  of  the  gradient.  In  §  4  we 
extend  the  method  to  multivariable  discrete  systems.  In  §  5  we  illustrate  the 
performance  of  the  algorithm  with  some  examples. 


2.  Computation  of  the  gradient  vector 

In  this  section  we  derive  explicit  expressions  for  the  cost  function  V  and  its 
gradient  vector  with  respect  to  the  coefficients  of  the  polynomials  a(z)  and  6(z). 
We  first  express  the  cost  function  in  terms  of  three  matrices,  each  of  which 
satisfying  a  certain  matrix  Lyapunov  equation.  These  Lyapunov  equations 
are  shown  to  admit  closed-form  solutions,  involving  differences  of  products  of 
triangular  Toeplitz  matrices.  Then  we  use  these  expressions  to  derive  a 
formula  for  the  gradient  vector. 


2.1.  The  cost  function 

Let  e(z)  be  the  z-transform  of  the  error  between  the  impulse  response  of  the 
given  transfer  function  and  that  of  the  reduced-order  approximate  model,  i.e. 


e(z) 


6(*) 

o(z) 


(7) 


Let  {Aj°,  0  <  i  <  qo  }  and  {h(,  0  ^  i  <  oo}  be  the  impulse-response  sequences  of 
I/«°(z)  and  l/a(z)  respectively,  i.e. 


1 


00 


00 


I  V“‘ 


(8) 


Using  these  sequences,  we  can  write  (7)  in  a  matrix  form.  We  shall  use 
semi-open  brackets  to  denote  semi-infinite  vectors  and  matrices,  so  that 


V 

'  V 

rvi 

^0 

ptl 

<* 

** 

V  V 

W  A.vV-'V 

— 

Ax  h0 

A„_,  A„_,.  .  .kq 

A. 

(9) 
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or  more  compactly  as 

(10) 

where  H°  has  dimensions  oo  x  A'  and  H  is  oox».  Assuming  that  both  a°(z) 
and  a(z)  are  stable,  all  semi-infinite  entities  in  (10)  have  finite  norms.  There¬ 
fore,  we  can  express  the  squared  norm  of  •  as 

F(a,  b)  A  *Te  =  bT Rb  -  2b0TQb  +  b^Sb0  (11) 

where 

RAITH-,  QAH^H;  SAH<”Ho  (12) 


The  dimensions  of  R,  Q  and  S  are  n  x  n,  X  x  n  and  X  x  X  respectively.  Let 
p(z)  be  a  monic  polynomial  of  degree  m,  where 

pto^l+P!*-1-!-  ...  +pmz~m 

We  define  C(p)  as  the  companion  matrix  of  p(2),  i.e. 

\~P\  -Pv-Pm 

0 

(13) 

1 

1 

The  matrices  R,  Q  and  S  can  be  characterized  in  terms  of  the  companion 
matrices  of  a(z )  and  a°(z )  as  follows. 

Lemma  1 

Each  of  the  matrices  R,  Q  and  S  is  the  unique  solution  of  a  matrix  Lyapunov 
equation 


R-C{a)RCT(a)  =  Enxn 

(14  a) 

Q-C(a°)QC*{a)~E„XH 

(14  ft) 

S-C(a°)SCr(a°)  =  ExxX 

(14c) 

Ekxl  is  a  matrix  of  dimension  lex l  haring  1  in  its  (1,  l)th  entry  and  zeros 
elsewhere. 

The  proof  is  by  a  direct  substitution  using  the  defining  relationships  (8). 
Existence  and  uniqueness  of  the  solutions  are  guaranteed  by  the  stability  of 
n(z)  and  a°(z)  (Lancaster  1969). 

The  next  lemma  gives  an  explicit  expression  for  the  solution  to  a  matrix 
Lyapunov  equation  of  the  type  appearing  in  Lemma  1. 

Lemma  2 

Let  p(z)  and  q(z)  be  two  stable  monic  polynomials  of  degree  m,  and  let  A' 
satisfy  the  matrix  Lyapunov  equation 


1 

C(p)  = 

0 


X  -  C(p)XCT(q)  = 

^mxm 


(15) 
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Then  : 


(1)  X  is  a  Toeplitz  matrix  (in  general,  non-symmetric) ; 

(2)  X~l  is  equal  to  QtP,r  -  PiQtT<  defined  as 


'  1 

'1  Pl  -Pm-l' 

o 

; 

ft  1 

•  Pi 

I  *,  *, 

o 

.fm-I-  fl  1. 

i 

‘  Pm 

9m  9m— !••• 

o 

Pm—l  Pm 

^in  tfm—l 

:  \ 

o 

_  Pi  ‘"Pm— l  Pm„ 

(16) 


The  lemma  can  be  proved  by  rather  tedious  algebraic  manipulations  or, 
more  easily,  by  using  results  from  the  theory  of  bi-orthogonal  polynomials  on  the 
unit  circle — see,  for  example,  Kailath  et  al.  (1978)  for  a  detailed  disoussion. 

To  use  Lemma  2  for  eqn.  (15  6),  a  slight  modification  of  this  equation  is 
necessary,  since  n  =deg  a  <  deg  o°  =  A'.  We  redefine  a(z)  as 

a(z)  =  1  +a12-1+  ...  +anz-"  +  0z-B_1  +  ...  +0z~‘v  (17) 

The  polynomial  d(z)  is  only  formally  different  from  a(z),  i.e.  d(z)=a(z)  for  all 
numerical  values  of  z.  Let  Q  be  the  N  x  N  matrix  satisfying  the  matrix 
Lyapunov  equation 

<3-C(a<»)<3(7T(o)  =  £VxA.  (18) 

Then  it  can  be  verified  that  Q  and  Q  coincide  in  their  n  leftmost  columns,  or  in 
other  word.*  <3  is  a  Toeplitz  extension  of  Q  to  a  square  matrix. 


2.2.  The  gradient 

We  now  have  all  the  necessary  relationships  for  computing  r  V fra  and 
dV/cb.  The  computation  proceeds  as  follows. 

ii  =  bTJ?b-2bOT^b  =  bT^b-2b«fr^6  (19) 

oa{  daf  dai  oai 

where  £  is  an  extension  of  b  to  dimension  S  by  adding  S  —  n  zeros.  Next  recall 
the  following  relationships  between  derivatives  of  matrices  and  of  their  inverses 


oaf  oai  cai  cai 

We  now  use  Lemma  2  to  express  R~l  and  Q~l  as 

Q~l  =  /l,®  A,r 


(20) 


(21  a) 
(21  6) 


Output-error  model  reduction 


101 


where  .4,,  .4,,  .i„  .J2,  Aj°  and  At°  are  defined  in  an  obvious  manner.  Differen¬ 
tiating  with  respect  to  a,  we  get 


3  R-^cA,  3^1'r_3£  3 A£_ 

3 a,  cai  1  1  3a,  3a,-  4  *  cai 


dQ~l  dA1 


=  -pA<"-At 


ca oa,- 


dAtT 

3a, 


(22  a) 

(22  6) 


Substituting  (20)  and  (22)  into  (21),  we  get 
0fzt  \  cai  da{  ) 


+  2bOT0  — 1  A.0  QB  -  26t<9t  ^  Ator  QT b“  (23) 

cai  ca, 


3X 


It  is  convenient  to  introduce  the  following  vectors 

r  —  Qrb° ;  s  =  QB  =  Qb  ;  t  =  Rb 
v1=»,41T.ftb;  v2  =  /liTi2b 
w,  =  At0  Qh  ;  w4  =  At°  QTb° 


and  then 


iI«2(-tT±iiy1  +  tTM*v2+rT^!w1-sT^Wt} 

cai  (  oa,  3a,  ca,  3a,  j 


(24) 


(25) 


Finally  we  need  an  expression  for  the  derivatives  of  the  matrices  ^4,,  .£„  At 
and  .4j.  Define  an  m  x  m  matrix  Z,nk  by 


(Zm% 


f1  =  »W-* 

I  0  ;  otherwise 


(26) 


Then  it  can  be  easily  verified  that 


3.4, 
3  a, 


dAt 

3a, 


=  Z 


,v 


oA,  _  ..  . 

_ ‘  _  7  .V-» 

3a, 


(27) 


Using  these  expressions  in  (25)  for  1  $  i  $  n  and  stacking  the  results  in  a  column 
vector  gives  the  desired  expression  for  3  l'/3a  as 


•••  l„_i  t„  o 

r 

' 

^n-l 

tn-l  \ 

o 

‘n  '  / 

O 

v,  +  2 

*2 

_  0 

.  <1  . 

. *»-. 

°  v  . 
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r 

*X 

.*  .* 

o 

4  * 

+  2 

/ 

w,-2 

:  \ 

W* 

/  /  O 

!  \  \ 

_rn+V-'rN  0 

*N-n+VaN-l  *X  0...0 

Finally.  dV/db  can  easily  be  computed  to  be 


dv 

rb 


=  2Rb  —  24?Tb°  =  2t  —  2 


(29) 


Lr»J 

In  the  next  section  we  show  how  the  computation  can  be  implemented  in  an 
efficient  manner  and  describe  other  components  of  the  output  error  algorithm. 


3.  Implementation  of  the  method 

In  this  section  we  discuss  several  issues  pertaining  to  the  implementation 
of  the  proposed  model-reduction  procedure.  First  we  consider  the  initialization 
problem  and  show  how  to  obtain  a  stable  reduced-order  initial  denominator 
polynomial.  Then  we  discuss  a  fast  computational  procedure  for  the  gradient 
vector.  Next  we  consider  the  problem  of  stability  monitoring  and  finally 
describe  the  gradient  search  procedure. 


3.1.  Initialization 

The  error  surface  corresponding  to  the  output-error  rational  approximation 
method  has,  in  general,  several  local  minima.  As  any  gradient  method  is  only 
guaranteed  to  converge  to  a  local  minimum,  it  is  imperative  to  choose  an  initial 
condition  which  is  sufficiently  close  to  the  global  minimum.  Furthermore,  it  is 
necessary  to  choose  a  stable  initial  condition  for  a(z)  and  to  keep  monitoring  the 
stability  as  the  search  proceeds. 

It  has  been  suggested  in  the  past  to  use  an  equation-error  approximation  of 
the  given  model  as  an  initial  condition  (Sanathanan  and  Koemer  1963). 
Unfortunately,  equation-error  approximations  are  not  guaranteed  to  be  stable. 
A  trivial  stable  initial  condition  is  a(z)  =  I,  but  this  may  be  too  far  from  the 
global  minimum  to  guarantee  convergence  to  this  minimum. 

We  propose  choosing  the  initial  a(z)  as  the  nth-order  orthogonal  polynomial 
of  the  given  a°(z)  on  the  unit  circle  (Szego  1967).  In  other  words,  a(z)  is 
defined  as  the  unique  solution  of  the  normal  equation 


‘  I  ‘ 

I 

** 
- 1 

_  6  _ 

(30) 


where  S„„l  is  the  (n+l)x(w  +  l)  principal  minor  of  the  Toeplitz  matrix  S 
defined  in  (12).  The  polynomial  a{z)  thus  defined  has  the  following  properties 
(Szegd  1967) 
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( 1 )  it  is  guaranteed  to  be  stable  ; 

(2)  it  is  an  optimal  nth-order  approximation  to  a°(z)  in  the  sense  that  it 
minimizes  the  prediction  error 

±1  I  a  (exp  ( jc»))  *  ^ 

2ni.\  a0  (exp  (jut))  w 

(3)  it  can  be  efficiently  computed  using  I-evinson’s  algorithm,  requiring 
about  n*  operations  (Kaiiath  1974). 

These  properties  make  a(z)  given  by  (30)  an  especially  attractive  choice  for  an 
initial  condition,  even  though  there  is  no  guarantee  this  will  lead  to  a  global 
minimum. 

As  we  shall  now  describe,  the  polynomial  h(z)  is  determined  at  each  iteration 
by  forcing  oVjcb  to  zero.  Hence  no  initial  condition  for  b(z)  is  required  here. 


3.2.  Efficient  computation  of  the  gradient 

As  we  have  shown,  both  the  cost  function  and  the  gradient  vector  require  the 
solution  of  matrix  Lyapunov  equations  of  the  form  of  (15).  Specifically,  (14  c) 
needs  to  be  solved  only  once,  while  (14  a)  and  (18)  need  to  be  solved  at  each 
iteration.  Consider  (15)  and  the  explicit  formula  (16)  for  the  inverse  of  its 
solution.  The  matrix  X  can  be  obtained  by  a  direct  inversion  of  the  right-hand 
side  of  (16),  but  this  would  require  m*  operations.  A  more  efficient  method  for 
inverting  this  matrix  is  by  the  so-called  inverse  Levinson  algorithm.  This 
algorithm  computes  the  UDL  decomposition  of  X~l  in  about  2m*  operations  (an 
operation  is  defined  here  as  one  multiplication  and  one  addition).  As  X  is 
Toeplitz,  it  is  fully  determined  by  its  first  and  last  columns,  and  those  can  be 
readily  computed  from  the  L-D-U  factors  of  X~l.  We  give  below  a  summary 
of  the  inverse  Levinson  algorithm,  skipping  the  proof  (see,  for  example,  Vieira 
and  Kaiiath  (1977)  for  the  symmetric  case). 


Step  1 .  Set 


?«,.=?<>  0<Um;  dm=  1 


(31) 


Step  2.  For  i  *■  m  down  to  i «  1,  do 

1 

Pi  =  -Put’  t<  =  T - — 

1  -  PfOf 


0  Pi-M-I  •••  Pi- t,l  1 

.*  ?<- i.i  0. 


"1 

Pi' 

Pi,  i  Pi,  i— l 

1 

£ 

l 

.  1  <7i,  i 

•  ?i.i- 1  ?i.ij 

(32  «) 


(32  h) 


(32  c) 
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Step  3.  Solve  the  following  equations  for  x,  and  x„,  the  first  and  last  columns 
of  X  respectively 


1 

r<q 

Pui  l'  O 

0 

•  • . 

xt  = 

_o_ 

*  1m- 1,1  ■ 

•  Im—Um—l 

X.  =* 

'o' 

1 

o 

lui 

1 

ft 

6 

do. 

(33) 


In  counting  the  number  of  operations  for  solving  the  two  Lyapunov 
equations  (14  a)  and  (18),  we  note  the  following  : 


(1)  Equation  (14  a)  is  symmetric  ;  hence pitj  =qu,  p(  =<r(  and  x,  is  sufficient 
to  determine  X.  Thus  the  total  count  of  operation  for  this  equation  is 
l-5n*. 

(2)  Equation  (18),  while  non-symmetric,  is  ‘  sparse  ’  in  the  sense  that  the 
corresponding  q(z)  polynomial  is  of  degree  n,  rather  than  AT.  The  total 
count  of  operations  taking  advantage  of  this  sparseness  is  about 
1-5AT*  +  l-5n*. 


The  gradient  search  procedure  can  be  improved  by  forcing  the  component 
dVjd b  to  zero  at  each  iteration.  This  has  the  effect  of  conditionally  optimizing 
V  with  respect  to  b  at  each  iteration  (where  tne  conditioning  is  on  the  current 
value  of  a),  thus  reducing  the  number  of  free  parameters  from  2 n  to  n.  As  we 
see  from  (29),  this  achieved  by  setting  b  to 


b=R-lQTbO  =  I{-l 


(34) 


The  computation  of  r  takes  X*  operations  and  the  solution  of  (34)  takes  2n* 
operations  (e.g.  by  substituting  for  R~l  its  expression  given  in  (16)).  The 
computation  of  t  is  then  saved,  since  now  t  =  [r, ...  r„]T. 

The  total  count  of  operations  can  now  be  computed  to  be  about 
4-5Ar*  +  2iVn  +  6nl.  This  can  certainly  be  considered  as  efficient;  by  com¬ 
parison,  a  more  conventional  solution  (say  of  the  form  used  by  Yahagi  (1981)) 
would  require  a  number  of  operations  proportional  to  nA'*. 


3.3.  Stability  monitoring 

Stability  monitoring  can  be  done,  in  principle,  by  solving  for  the  roots  of 
a(z)  and  checking  that  they  are  all  inside  the  unit  circle.  This,  however,  is  an 
undesirable  approach,  since  it  significantly  increases  the  computational  burden 
if  the  degree  of  o(z)  is  relatively  large.  Alternatively,  the  stability  of  a(z)  can 


V 


>/■ 
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be  tested  by  the  Schur-Cohn  test  (Jury  1974),  which  does  not  require  a  factori¬ 
zation  of  the  polynomial.  An  interesting  feature  of  our  algorithm  is  that  it,  in 
fact,  includes  a  stability  test.  The  solution  of  (14  a)  via  the  algorithm  (31)-(33) 
is  equivalent  to  the  Schur-Cohn  test  in  the  symmetric  case  (Vieira  and  Kailath 
1977).  The  condition 

|p,|<l,  t«l,...,n 

is  necessary  and  sufficient  for  stability  of  the  given  polynomial.  Thus  our 
algorithm  provides  stability  monitoring  at  no  extra  computational  cost. 

3.4.  The  gradient  search  procedure 

Once  a  closed-form  expression  for  the  gradient  is  available,  one  of  many 
existing  gradient  methods  can  be  used  for  minimizing  the  error  function  K(a,  b). 
We  have  chosen  to  use  the  Fletcher-Powell  method  (Luenberger  1973),  known 
for  its  excellent  convergence  rate  and  relative  ease  of  implementation.  An 
important  part  of  this  method  (as  well  as  of  virtually  all  gradient  methods)  is 
the  line  search  procedure,  namely,  a  search  for  a  local  minimum  of  the  cost 
function  at  the  direction  used  at  each  iteration,  as  a  function  of  the  step  size. 
In  our  case,  a  certain  difficulty  occurs  due  to  the  fact  that  the  a  vector  is 
constrained  to  be  in  the  open  set  Q  =  {a  :  a(z)  is  stable}.  On  the  boundary  of 
this  set  the  error  V(a,  b)  approaches  infinity,  and  it  is  not  defined  outside  the 
set.  Thus,  whenever  the  poles  of  a°(z)  are  near  the  unit  circle,  great  care  is 
needed  in  performing  the  line  search  to  stay  within  the  permitted  region  Q. 
We  have  found  the  golden  section  search  procedure  (Luenberger  1973)  very 
useful  in  this  case,  since  it  uses  the  values  of  the  cost  function  only  for  magnitude 
comparisons  and  makes  no  use  of  derivatives.  Thus,  by  assigning  very  high 
cost  to  an  unstable  a  (say  near  the  value  of  the  computer  overflow)  the  line 
search  can  be  forced  to  yield  only  stable  values  of  a. 


4.  Model  reduction  of  multivariable  systems 

In  this  section  we  consider  the  case  where  both  g°{z)  and  g(z)  are  p  x  m 
transfer-function  matrices,  rather  than  scalars.  A  natural  rational  description 
of  such  matrices  is  in  terms  of  so-called  matrix  fraction  descriptions  (MFD) 
(Kailath  1980).  Formulating  the  output-error  model-reduction  problem  in 
terms  of  MFDs  leads  to  Lyapunov  equations  in  block-companion  form.  Unfor¬ 
tunately,  such  equations  do  not  appear  to  admit  closed-form  solutions  of  the 
type  shown  in  (16).  (It  is  worth  noting  that  substituting  matrices  for  scalars 
in  ( 16)  does  not  lead  to  a  correct  solution  of  (15)  in  the  matrix  case.)  Therefore 
we  have  chosen  not  to  use  MFD  representations  here  but  take  a  different 
approach. 

Let  a°(z)  be  the  characteristic  polynomial  of  the  system  whose  transfer- 
function  matrix  is  g°(z).  Then  g°(z)  can  be  written  as 


<Az) 


B°(z) 

a°(z) 


V  z~l+  ...  +Bn°z-* 

1  -t-®,0  z~l  +  ...  +aN°z~x 


(35) 


where  {Bj0,  ...,  fly0}  are  px  m  matrices.  As  before,  g°(z)  is  assumed  to  be 
stable  and  strictly  proper. 
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We  wish  to  approximate  g°(z)  by  the  nth-order  pxm  transfer-function 
matrix 


B(z)  Bxz~l+  ...  +  Bnz~n 
a(z)  l+a,2~l+  ...  +anz~n 


(36) 


We  proceed  as  in  §  2.1,  expressing  the  error  as  a  semi-infinite  vector  and 
then  expressing  F  as  the  lt  norm  of  this  vector.  It  will  be  convenient  to 
introduce  the  following  notation  :  let  B(°  be  a  row  vector  of  dimension  pm, 
obtained  from  Bt°  by  stacking  its  rows  in  their  natural  order.  B{  is  defined  in  a 
similar  manner.  Using  these  definitions,  we  can  express  the  error  vector  as 


r  a»°  l 

®1° 

*  K  o' 

bi 

:  \0 

; 

«* 

= 

* 

— 

• 

A.v-10 . -V 

i . Aq 

or  more  compactly  as 

e  =  H°B°  -  HB 


(37) 


(38) 


The  element  et  of  e  is  now  a  row  vector  of  dimension  pm  whose  entries  are  the 
pm  components  of  the  impulse  response  at  time  ».  The  cost  function  F  is  now 
given  by 

F(a,  B)4tr  {®T*}=tr  {BT5B  — 2B«TQB  +  B«TSB0}  (39) 

where  tr  {  - }  denotes  the  trace  operator  and  R,  Q  and  S  are  defined  as  in  §  2.1. 

The  rest  of  the  procedure  is  similar  to  the  one  given  in  §  2,  with  some  minor 
modifications.  In  particular : 

( 1 )  the  matrices  R,  Q  and  8  are  obtained  exactly  as  before  ; 

(2)  the  gradient  component  dV/oB  is  now  given  by 

dV/dB  =  2RB-2QTB°  (40) 

By  setting 

B  =  /?-’QtB°  (41) 

at  each  iteration,  the  dimensionality  of  the  problem  decreases  from 
{mp  +  l)w  to  n.  This  entails  a  considerable  saving  in  the  amount  of 
operations  ;  therefore  it  is  highly  recommended  here. 

(3)  The  expression  (25)  for  SF/Sa,  basically  remains  the  same,  except  for 
the  need  to  take  the  trace  of  the  right-hand  side.  This  causes  some 
difficulty  in  obtaining  an  expression  of  the  form  (28)  for  3F/3 a.  How¬ 
ever,  such  an  expression  is  not  really  needed  for  practical  implementation 
of  the  method.  It  is  sufficient  to  compute  the  dV/dat  individually, 
and  then  stack  them  in  a  column  vector  of  dimension  n. 


5.  Numerical  results 

In  this  section  we  demonstrate  the  performance  of  the  algorithm  by  some 
numerical  examples.  We  have  chosen  to  test  transfer  functions  which  have 
poles  near  the  unit  circle,  as  these  cases  are  usually  quite  difficult  to  handle. 
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The  first  case  uses  the  8th-order  transfer  function  whose  denominator  and 
numerator  polynomials  are 

a°(z)  -  1  -  4-082z‘l  +  7-226&S-*  -  6-4408Z"*  +  1  SlOSz-4 

+  2-0443z“*  -  2-4197z-«  +  1  0356z-7  -  0-  I516z~*  (42  a) 

b°(z)  «  z-‘  -  0-357Z-*  +  0-2038Z-*  -  0-0848z~*  -  0-0493 z~*  -  0-1 92z~*  (42  b) 

The  poles  of  this  system  are  at  -0-7,  0-9274  ±  jO-30 15,  0-7883  ±  j0-5730,  0-3, 
and  0-3827  ±  j0-7867. 

The  order  taken  for  the  reduced  model  was  n  =  4.  Figure  2  shows  the 
impulse  response  of  the  approximation  corresponding  to  the  initial  choice  of 
a(z)  as  described  in  §  3. 1 ,  against  the  impulse  response  of  the  full  model.  Figure 
3  shows  the  corresponding  frequency  responses.  The  poles  of  the  initial 
reduced-order  model  are  at  0-9205  ±  j0-3213  and  0-7754  ±  j'0-6058.  We  see 


Example  I  —  Impulse  response  of  initial  approximation 


REDUCED 


Figure  3.  Example  1 — Frequency  response  of  initial  approximation. 
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that,  even  though  the  poles  are  rather  close  to  the  dominant  poles  of  the  full- 
order  model,  this  approximation  is  nevertheless  poor.  The  squared  error  of  this 
approximation  is  728-6,  or  about  25%  of  the  squared  impulse  response  of  the 
full  model,  which  is  2865*2.  Figures  4  and  5  show  the  impulse  and  frequency 
responses  respectively  of  the  approximation  obtained  after  eight  iterations  of 
the  algorithm.  The  transfer  function  of  this  approximation  iB 

b(z)  —  0-3945z-1  +  4-  83 1 82-8  —  4-8 1 69z-3  4-  0-9268Z-4 

o(2)  "  1  -  3-4301z-»  +  4-8228Z-*  -  3-25992"*  +  0-903 lz~* 

The  poles  are  at  0-9271  ±j'0-3099  and  0-7879  ±j'0-5741.  The  squared  error  is 
8-48,  i.e.  about  1%  of  its  initial  value  !  The  match  of  the  impulse  responses  is 
excellent.  The  match  of  the  frequency  responses  is  excellent  down  to  20  dB 
and  then  starts  deteriorating.  This  is  an  obvious  result  of  the  fact  that  the 
output-error  method  weighs  the  error  linearly,  while  the  frequency  response  is 


9.  49.  39.  129.  169  .  299. 


- -  REDUCED 

Figure  4.  Example  1  —  Impulse  response  of  final  approximation. 


REDUCED 


Figure  5.  Example  1 — Frequency  response  of  final  approximation. 


Output-error  model  reduction 


109 


shown  on  a  logarithmic  scale.  Thus  one  should  not  expect  a  good  match  of  the 
frequency  responses  at  frequencies  where  the  energy  density  is  low. 

The  second  example  uses  the  lOth-order  transfer  function  when  denominator 
and  numerator  polynomials  are 

o°(z)  -  1  -  4-4I58s-i  +  8-4582z“*  -  8-2398z“4 

+  2-5658Z-*  +  3-2817Z-*  -  4-2475z~«  +  l-5203z-» 

+  0-681 97z~*  -  0-83 1 62z~»  +  0-25 1 5  lz~ 10  (44  a) 

h°(s )  -  -  0-31272z-‘  -0-39268Z-*  +  2-3363z~* 

-  2-031 8z~<  +  0-807632 -  0-6027z~«  -  0-86225z-J 

+  1  -6256z~*  +  0-03 1 42z~*  —  0-64564z-10  (44  6) 

This  example  is  taken  from  Kung  and  Lin  (1980).  The  poles  of  this  system  are 
at  0-9561  ± j0-2721,  0-3827  +  j0-7867,  -  0-6349  ±j0- 1454,  0-871 1  ±  jO-4517, 
and  0-6329  +  JO-6430.  The  reduced  order  was  taken  to  be  n  =  6. 


- -  REDUCED 

Figure  6.  Example  2 — Impulse  response  of  initial  approximation. 


-  REDUCED 

Figure  7.  Example  2 — Frequency  response  of  initial  approximation. 
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Figures  6  and  7  show  the  impulse  and  frequency  responses  of  the  initial 
approximation.  The  squared  error  of  this  approximation  is  295-5,  while  the 
energy  of  the  full  model  is  778-9.  Again,  the  initial  approximation  is  definitely 
poor  in  this  case. 

Iterating  the  algorithm  22  times  gave  rise  to  the  approximation 
-  0- 1 1 12-‘  -  1  -4974s-*  +  3-85782“* 

b(z)  —  4-32482-* +  4-9636Z'4  —  3-01 19s"* 

a(z)  1  -  3-74732-*  +  6-68452“*  -  7-26822-*  '  ' 

+  5-09342-*  -  2-21 1 52-*  +  0-4930z-« 

The  poles  are  at  0-9560  ±>0-2722,  0-5922  ±  j0-639l,  0-3254  ±j0- 7425.  Some¬ 
what  surprisingly,  four  of  the  poles  are  not  very  close  to  the  dominant  poles  of 
the  full-order  system.  The  approximation  is  still  very  good,  as  is  shown  in 
Figs.  8  and  9.  The  squared  error  of  the  approximation  is  23-7,  or  about  8%  of 
its  initial  value. 


REDUCED 


Figure  8.  Example  2 — Impulse  response  of  final  approximation. 


-  REDUCED 

Figure  9.  Example  2 — Frequency  response  of  final  approximation. 
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The  last  example  again  uses  the  model  (44)  but  takes  the  reduced  order  as 
n  =  4.  The  initial  approximation  is  shown  in  Figs.  10  and  11,  and  the  final 
approximation  in  Figs.  12  and  13.  The  initial  and  final  squared  errors  are 
498-9  and  52-5  respectively.  The  number  of  iterations  needed  to  reach 
convergence  was  eleven.  The  resulting  approximate  model  is 

6(z)  2-64892-*-  10-96802-*+ 15-52052-*- 7-2175Z-* 

a(2) =  1  -  3-18562-*  +  4-26162-*  -  2-86 1 52“*  +  0-82852"4  (46) 

The  poles  are  at  0-9561  ±  j0-2724  and  0-6367  ±  jO-6579. 

The  reader  is  referred  to  Kung  and  Lin  (1980)  for  a  comparison  with  other 
approximation  methods  (singular-value  decomposition,  Hankel-norm  approxi¬ 
mation,  and  dominant-mode  approximation). 


-  REDUCED 

Figure  10.  Example  3 — Impulse  response  of  initial  approximation. 


-  REDUCED 

Figure  11.  Example  3 — Frequency  response  of  initial  approximation. 
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I.  44.  84.  12*.  16*.  a**. 


4  .  44  .  84.  128.  164.  288. 


-  FULL 

-  REDUCED 

Figure  12.  Example  3 — Impulse  response  of  final  approximation 


8. 888  8.188  8.288  8.388  8.488  8.588 


8.888  8.188  8.288  8.388  8.488  8.588 


-  FULL 

-  REDUCED 

Figure  13.  Example  3 — Frequency  response  of  final  approximation. 


6.  Conclusion 

We  have  presented  an  efficient  algorithm  for  output-error  model  reduction 
of  linear  discrete  systems.  The  key  step  was  the  development  of  a  closed-form 
expression  for  the  gradient.  Furthermore,  this  expression  was  shown  to  be 
efficiently  computable,  requiring  a  number  of  operations  proportional  to  »VS, 
where  X  is  the  order  of  the  given  model.  The  proposed  optimization  technique 
employs  a  special  initialization  procedure  and  the  stability  of  the  algorithm  is 
guaranteed  by  a  built-in  stability  test. 

The  method  can  easily  be  extended  to  model  reduction  of  multi-input- 
multi-output  systems.  Somewhat  surprisingly,  this  is  achieved  without  the  use 
of  matrix  fraction  descriptions,  by  dealing  with  the  characteristic  polynomials 
directly. 

The  new  algorithm  has  potential  applications  to  filter  design,  system 
identification  and  control  systems.  Our  particular  motivation  in  developing 


Output-error  model  reduction 


113 


the  algorithm  was  for  possible  application  to  adaptive  control  of  large-scale 

systems.  Some  results  obtained  in  this  specific  application  are  reported  in 

Friedlander  and  Porat  (1982). 

References 

Aplevich,  J.  D.,  1973,  Int.  J.  Control,  17,  565. 

Deczky,  A.  G..  1972, I.E.E.E.  Trans.  Audio  Elect roacoust.,  20,  257. 

Friedlander,  B.,  and  Porat,  B.,  1982,  Proc.  21s<  I.E.E.E.  Conf.  Decision  and 
Control,  Florida,  U.S.A. 

J  cry,  E.  S.,  1974,  Inners  and  Stability  of  Dynamic  Systems  (Chichester  :  John  Wiley). 

Kailath,  T.,  1980,  Linear  Systems  (Englewood  Cliffs,  N.J. :  Prentice  Hall) ;  1974, 
I.E.E.E.  Trans.  Inf.  Theory,  20,  146. 

Kailath,  T„  Vieira.  A.,  and  Morf,  M.,  1978,  SIAM  Rev.,  20,  106. 

Kuno,  S.,  1980,  Proc.  Joint  Conf.  Automatic  Control,  Paper  No.  FA8-D,  San  Francisco, 
CA,  U.S.A. 

Kuno,  S.,  and  Lin,  D.  \V.,  1980,  Reduced  Order  System  Modeling  via  Singular  Value 
Analysis  (University  of  Southern  California,  Dept.  Electrical  Engineering). 

Lancaster,  P.,  1969,  Theory  of  Matrices  (New  York  :  Academic  Press). 

Landau,  Y.  D.,  1979,  Adaptive  Control :  the  Model  Reference  Approach  (New  York  : 
Marcel  Dekker). 

Luenberger,  D.,  1973,  Introduction  to  Linear  and  Nonlinear  Programming  (New 
York  :  Addison  Wesley). 

Moore,  B.  C.,  1978,  Proc.  I.E.E.E.  Conf.  Decision  and  Control,  New  Orleans,  LA, 
66-73. 

Sanathanan,  C.  K.,  and  Koernbr,  J.,  1963,  I.E.E.E.  Trans,  autom.  Control,  3,  56. 

Steiolitz,  K.,  and  McBride,  L.  E.,  1965,  I.E.E.E.  Trans,  autom.  Control,  10,  461. 

Szego,  G.,  1967,  Orthogonal  Polynomials  (Providence,  RI :  American  Mathematical 
Society),  Colloquium  Publications  No.  23,  3rd  edn. 

Vieira,  A.,  and  Kailath,  T.,  1977,  I.E.E.E.  Trans.  Circuits  Syst.,  24,  218. 

Yahaoi,  T.,  1981,  I.E.E.E.  Trans.  Acoust.,  Speech,  Sig.  Processing,  29,  245. 


v.vy 


APPENDIX  C 

AN  OUTPUT  ERROR  METHOD  FOR 
REDUCED  ORDER  CONTROLLER  DESIGN 


-.'TVr 


IEEE  TRANSACTIONS  ON  AUTOMATIC  CONTROL,  VOL  AC-29.  NO.  7.  JULY  1984 

Technical  Notes  and  Correspondence. 


An  Output  Error  Method  for  Reduced  Order 
Controller  Design 

BOAZ  PORAT  and  BENJAMIN  FRIED  LANDER 

Abstract — An  efficient  computational  tcchmquc  is  presented  for  the 
design  of  reduced  older  controllers  for  Knew  discrete-time  systems.  The 
technique  is  based  on  the  minimi?  artna  of  the  output  error  between  the 
dosed-toop  system  and  a  specified  reference  model. 

I.  Introduction 

This  note  is  concerned  with  the  problem  of  designing  reduced  order 
controllers  for  discrete  control  systems  using  a  least  squares  error  crite¬ 
rion  with  respect  to  a  given  reference  model.  Let  the  plant  under  consider¬ 
ation  be  represented  by  the  transfer  function 

C(-)-/(i)/*(0“  (  i>*")  j Ls.-*v")  :  (1) 

also  let  the  reference  model  be  represented  by  the  transfer  function 

G*(;)-w(i)u(r)-  (  L  ’V*'1)  / I  “fi"*’)  •  (2) 


Both  the  plant  and  the  reference  model  are  assumed  .0  be  strictly  proper. 
The  reference  model  is  assumed  to  be  asymptotically  stable.  Denote  by 
H(:)  the  transfer  function  of  the  desired  nth  order  cascade  compensator, 
where 

//(.-)-*(. •)/<.(.*)-(  £  *ia— J  j  £«,--*-)•  (3) 

The  closed-loop  transfer  function  <7r( z )  is  given  by 

G,(z)-t>(z)f(z)/(a(z)g(z)  +  b(z)/(z))*c(z)/d(z).  (4) 

Let  <(:)  denote  the  difference  between  the  transfer  functions  of  the 
reference  model  and  the  actual  closed-loop  system,  i.e.,  <(;)-(/*<;)- 
C,  ( : ).  The  cost  function  to  be  minimized  is  given  by 

K(<i1.  --.a..b0.  --.b.)-^/jr(e^) \ldu>-  £  (5) 


The  aim  is  to  find  a  cascade  compensator  h(z)/a{z)  that  will  minimize 
V.  It  will  be  assumed  that  n  <  N  +  Af,  so  that  V  -  0  is  impossible  in 
general.  Previous  work  on  this  problem  includes  (l|-(3].  The  contribution 
of  the  present  work  is  an  efficient  computational  scheme  for  the  gradient 
vector  of  V  with  respect  to  the  coefficients  of  a(z)  and  !>(:).  U sing  some 
facts  from  theory  of  discrete  Lyapunov  equations  and  Toeplitz  matrices, 
we  derive  an  algorithm  for  computing  the  gradient  vector  in  a  number  of 
operations  proportional  to  (Af  +  n)2.  The  gradient  is  then  used  in  a 
Fletcher- Powell  minimization  algorithm  to  optimize  the  controller 
parameters. 
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IL  Computation  of  the  Cost  Function  and  the  Gradient 

Recall  that  the  error  between  the  transfer  functions  of  the  reference 
model  and  the  closed-loop  system  is  given  by 

<(z)~"(z)/u(z)-c(z)/d(z).  (6) 

The  polynomials  c(z)  and  d(z)  are  of  degrees  a  +  jV  - 1  and  S  *  n, 
respectively,  while  w(r)  and  u(z)  are  of  degrees  M  - 1  and  M,  respec¬ 
tively.  It  is  convenient  to  multiply  both  w(r)  and  i<(z)  by  and 

redefine 


£  ulz!l*m~‘\  w(r)-  Z  "I*"”"-  (?) 

f-1  1-1 

Let  { A?:  0  <  /  <  oo }  and  { h,;  0  <  i  <  oo }  be  defined  by 
00  00 

zs*"/it(z)~  Z  *?*'*:  z"*Vd{z)~  Z  M"  (8) 

i—0  i—0 

Then  we  can  write  (6)  in  the  form 

(9) 

where  t  is  the  semi-infinite  vector  •  ••]r,  H°  and  H  are  lower 
trapezoidal  semi-infinite  Toeplitz  matrices  of  width  N  +  a  whose  first 
columns  are  ' lr  and  ■  •  •  ]r.  respectively,  *-(»■,■■• 

w...]r  and  c-  [c,  •  ■  •  cw.„]r.  Assuming  that  both  «(r)  and  dfz)  are 
Stable,  all  semi-infinite  entities  in  (9)  have  finite  norms.  Therefore,  we  can 
express  the  square  norm  of  <  as 

/  -  trt  -  ctHtHc  -  2wtHotHc  *  »TH0TH°w 
*•  cTRc  -  2wtQc  +  wTS»  (10) 

where  the  definition  of  R,Q,S  is  clear  from  (10).  Introduce  the  notation 
C(p)  for  the  top-row  companion  matrix  of  the  polynomial  p<z): 


C(F)- 


The  matrices  R,Q,S  can  be  expressed  in  terms  of  the  companion 
matrices  C(d)  and  C(u)  as  follows. 

Lemma  l:  Each  of  the  matrices  R,  Q,  and  5  is  the  unique  solution  of 
a  matrix  Lyapunov  equation: 

R-C(d)RCr(d)-£; 

Q  ~  C(u)QCr(  d )  •  £; 

5-C(u)5Cr(«)-£  (11) 

where  £  is  a  matrix  having  1  in  its  (1.1)  position  and  zeros  elsewhere.  The 
proof  is  by  a  direct  substitution  of  ti-e  definitions  of  R,  Q,  and  S  into 
(11).  Uniqueness  of  the  solution  is  guaranteed  by  the  stability  of  u< : )  and 
d(z).  An  algorithm  for  solving  these  equations  in  a  number  of  operations 
proportional  lo  (.V  +  a)2  is  given,  e.g.,  in  (4).  The  following  lemma  gives 
an  explicit  expression  for  the  solution  to  a  matrix  Lyapunov  equation  of 
the  type  appearing  in  Lemma  1. 

Lemma  .V  Lee  p(z)  and  qtz)  t/<  two  stable  polynomials  of  degree  m. 
and  let  X  satisfy  the  matrix  Lyapunov  equation 


C(p).VCr(q)-£. 
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Then  1)  X  is  a  Toeplitz  matnx  and  2)  AT 1  is  given  by 
X-'-QSl- hQl 

1  Ft  •  Pm- 1 
1 

0  p , 

1 

Pm  0 

Pm—\ 

P\  Pm - 1  Pm 


1  0 

?»- 1  •••  9t  1 


9i 

9m—  | 

9m  . 

03) 


O, +  +  (15) 

R“,"D|I>ir-  Oj2>f;  (?"* "  D|l/lr-lij2)f ;  S-'-UxU,r-UiU{ 

(16) 

We  now  have  all  the  necessary  relationships  lor  computing  ( 9V/9a„ 
1  <  i  <  « }  and  { 9v/9b„  0  <  i «  » ).  The  computatioa  proceeds  as  fol¬ 
lows; 


9o 

da, 

9V_ 

9b, 


eT^.c-lwT*2.c 

da,  da, 


Recall  that 


(17) 

(18) 


This  lemma  can  be  proven  using  results  from  the  theory  of  bionhogonal 
polynomials  on  the  unit-circle  (5).  It  will  be  used  here  to  write  down 
explicit  expressions  for  R~l.  Q~l,  and  S'1  as  follows.  For  a  polynomial 
p(c)  of  degree  +  let  th t(N  +  n)x(N  +  *)  matrices  Px  and  F2 
be  defined  by 


0  <  i  -  j  <  m 
otherwise 

otherwise. 


With  these  definitions  we  have 


a* 

9a, 


9£ 

da, 


(19) 


and  similarly  for  the  derivatives  with  respect  to  h,.  Hence 


£— Vg1**^*  tm 


To  get  explicit  expressions  for  the  derivatives  of  the  inverse,  let  us  define 
the  Ic-shift  matrix  of  dimensions  (N  *■  n)x(N  *  n)  by 
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Z*- 


{ 


1; 

0; 


» -  i  “  * 

otherwise. 


(») 


Thai!  is  easy  to 


that 


«2t 

9m, 


-r; 


ida 


-Z*- 


98, 

9b, 


(23) 


Using  (IS).  (16).  and  (23)  we  get 

^.-rc,j3,r+  0,cf(r)T- z— <ho[ - DfiHz—)r  (24) 

-  rFtf*  thF,r{2‘)T-  Z—FjD[-  DiF-[(Zt~')T  (23) 
^  -  Z?G\U\-  Ufih  Z-)  r  (26) 

Ql-rFtf-mFKz—)7  (27) 

Finally.  de/db,  is  ghieii  by 

r  (28) 

i  jy  a  -  i 


Equations  (20).  (21),  (24)-(2>)  fora  a  complete  algorithm  for  the  gradient 
vector.  It  can  be  checked  that  the  number  of  operations  involved  in  the 
computation  is  proportional  to  (N  +  a)1. 


Ill  Implementation  or  tub  aloomtkm 

The  algorithm  described  above  was  implemented  on  a  VAX  780.  in 
Fortran  77.  The  Lyapunov  equations  (11)  were  solved  using  a  nonsym- 
OKtric  version  of  the  inverse  Levinson  algorithm  (6}.  This  algorithm  was 
also  used  for  monitoriag  the  stability  of  die  d(z)  polynomial  at  each 
iteration.  The  optimization  method  used  was  Fletcher-Poweil  (7J.  with  a 
golden  section  line  search.  The  reason  for  this  particular  combination  was 
that  it  requires  only  one  gradient  computation  per  iteration,  thus  helping 
to  reduce  the  overall  number  of  computations. 


IV.  an  Example 

The  following  example  win  serve  to  illustrate  the  performance  of  the 
algorithm: 

G(z)  -  *4/(4*  -3-3r*  +4.68Z1  -2.921r:  +0.8263*  -0.083)  (29) 

C,(*)-0.2S*V(** -1-23* +0.5).  (30) 

The  order  a  of  the  cascade  compensator  was  chosen  as  2.  An  initial 
stability  compensator,  found  by  trial  and  error,  was  taken  as 

#(r)-(1.3s,-2r+0.823)/r2.  (31) 

Fig.  1  shows  the  impulse  respome  of  the  initial  closed-loop  system 
compared  to  that  of  the  reference  model  The  optimal  controller,  obtained 
after  23  iterations  of  the  algorithm,  was  found  to  be 

H(t)  -  (0.248rJ  -0.47r  +0.234)/(r:  +0116*  -0.033).  (12) 

Fig.  2  shows  the  impulse  response  of  the  final  dosed- loop  system  com¬ 
puted  to  that  of  the  reference  model.  In  this  example,  die  total  square 
error  of  the  optimal  solution  w as  about  0.03  percent  of  the  reference 
model  impulse  response  energy. 


controllen  was- presented.  The  algorithm  uses  inpul/output  (transfer 
function)  descriptions  of  the  plant,  the  reference  modd  and  the  controller, 
rather  than  the  mom  «— wiy  used  ime-apnee  descriptions.  The  two 
descriptions  arc  mathematically  equivalent,  but  lead  to  Afferent  computa¬ 
tional  procedures. 

The  relative  computational  efficiency  of  the  propoacd  algorithm  opens 
the  way  for  using  it  as  part  of  an  adaptive  reduced  order  controller  where 
the  procedure  needs  to  be  performed  repvauxfly.  This,  however,  is 
outside  the  scope  of  this  note. 
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Noopnnunetrk  Algorithm  for  Input  Signals 
Identification  in  Static  Distributed-Parameter 
Systems 

EWARYST  RAFAJLOWICZ 

Abuma— In  tMeemrijpandiari,  a  annpmmmrir  algoriduu  for  Idapti- 

propoaad  and  innatigmid.  Iwagral  mtan-aqnasu  cowurginf*  of  the  a%iK 
rithiu  la  proved  far  an  inflate  number  of  point  mmeanmentsef  the  system 

by  Rutkowdd  (t#i  far  iiuupwnmurtr  fanctian  fitting.  and  in  a  common 


I.  Introduction 

The  aim  of  this  correspondence  is  to  propose  and  investigate  an 
algorithm  for  identification  of  an  unknown  input  signal  or  an  excitation 
of  a  static,  linear  distributed-parameter  system  (DPS)  from  point  mea¬ 
surements  of  its  state.  Problems  of  this  type  arise  in  the  nr**  of  wafer  and 
air  pollution,  electromagnetic  beating,  vibrations  isolation,  etc.  and  have 
been  treated  by  several  authors  (mainly  from  a  computational  point  of 

view)l3H6V 

Theoretical  analysis  of  such  problems  is  a  difficult  task  since  they  are 
ill-posed  in  the  sense  pf  Hadantard  (2J.  [7J.  This  difficulty  is  usually 
avoided  by  assuming  a  priori  Uut  the  unknown  excitation  belongs  to  a 
certain  parametric  class  and  only  iu  parameters  are  estimated  [1|. 

In  (his  correspondence  no  assumptions  of  this  type  are  made,  and  thus, 
the  proposed  algorithm  is  a  nonparametric  one  Its  main  advantage  is 
asymptotic  optimality  (AO),  understood  as  the  integral  mean-square 
convergence  (IMS C)  as  the  number  of  measurements  approaches  to 
infinity. 

It  should  be  remarked  that  the  proposed  algorithm  is  a  modified 


V.  Conclusions 

A  computationally  efficient  algorithm  for  the  design  of  reduced  order 
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Lattice  implementation  of  some  recursive  parameter-estimation 
algorithms-* 

BENJAMIN  ERIEDLANDERJ 

Linear  dynamic  motielx  of  plant*  are  ueuaily  parametrised  by  tie*  i-uoffii-iout*  id 
difference  equations,  lattice  structures  ami  their  reflection  coefficients  provide  an 
alternative  paramet fixation  tliat  offers  several  advantage*,  including  numerical 
robustness,  computational  efficiency  anti  ease  of  liartlware  implementation.  Tlx- 
recursive  square-root  normalixetl  lattice  versions  of  the  following  well-known 
parameter-estimation  algorithms  are  presented  ;  recursive  least  -squares,  recursive 
instrumental  vartahle.  extemlml  least -squares,  ami  recursive  maximum-liki-hlsssl. 

1.  Inohctiw 

The  need  for  real-time  system  identification  led  to  the  development  of 
numerous  recursive  parameter-estimation  algorithms.  The  most  commonly 
used  algorithms  are  related  to  linear  input-output  models  described  by 
difference  equations  of  the  type 

.vs  xh 

j um-  £  a&i -<+  £  t  (i) 

l- 1  i-O 

where  «„  y,  denote  the  input  and  output  of  the  plant  and  t,  denotes  a  distur¬ 
bance  process.  Parametrizing  the  plant  by  the  coefficients  {a(,  bt]  of  the 
difference  equation  seems  to  be  a  natural  choice,  resulting  in  an  estimation 
problem  of  the  type  encountered  in  regression  analysis.  The  regression 

variables  in  this  case  are  simply  the  input  {u, . W/-.v«}  *nd  output 

{y,_, . y,_.v.i}  variables.  Other  parametrizations  are.  of  course,  possible. 

To  see  this,  rewrite  ( I )  in  the  form 

where 

-ift-i . ~ Ht-x A*  w . w,-.v/rlT 

lai . ns.v  ^i»’  —  fc.v»lT 

The  parameter  vector  8  and  the  vector  of  regression  variables  <f>,  can  he  replaced 
bv  3,*.Sd,  where  S  is  an  arbitrary  (possibly  time-varying)  non¬ 

singular  matrix,  since  clearly 

0  +  *,-$*  9  +  *,  (3) 

This  leads  to  an  infinite  number  of  possibilities  for  parametrizing  the  plant. 
An  interesting  choice  is  to  pick  8  so  that  the  regression  variables  will  he  un¬ 
correlated.  This  can  be  done  by  letting  be  the  lower-triangular  square-root 
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of  the  covarience  matrix  of  <f>t,  i.e. 

SST-E{M'T}  (4) 

in  this  case 

E{M?}=S-'E{M,T)S-r-l  (5) 

In  other  words,  the  original  set  of  regression  variables  is  replaced  by  a  Gram- 
Schmidt  orthogonalized  set.  Square-root  procedures  in  linear  least-squares 
estimation  are  known  to  have  good  numerical  properties  and  to  be  more  robust 
than  techniques  involving  the  covariance  matrix  itself  (Bierman  1977,  Lawson 
and  Hanson  1974).  Note  that  the  least-squares  estimate  of  the  parameter 
vector  S  is  given  by 

In  other  words,  the  parameters  can  be  interpreted  as  the  cross-correlation 
between  the  data  »/,  and  the  regression  variables.  Such  parameters  have  been 
used  for  quite  some  time  under  the  name  of  PARtial  CORrelation  (PARCOR) 
coefficients  in  the  analysis  of  time-series,  especially  in  speech  applications 
(Markel  and  Gray  1976).  This  parametrization  is  related  to  lattice  structures 
instead  of  the  tapped  delay -line  structure  inherent  in  the  difference  equation  (1 ). 

Lattice  forms  are  widely  used  in  signal-processing  applications  involving 
linear  filtering  and  prediction.  They  are  known  to  have  a  number  of  attractive 
features  including :  (i)  good  numerical  behaviour  on  finite-word-length 

processors ;  (ii)  an  orthogonality  (decoupling)  property  :  the  signals  pro¬ 
pagating  in  a  lattice  filter  are  uncorrelated.  (One  manifestation  of  this  property 
is  the  fact  that  when  the  filter  order  is  increased,  one  has  to  add  an  additional 
section  to  the  filter  without  changing  the  previous  sections.  In  other  words  the 
(*V  +  l)th-order  lattice  predictor  is  the  same  os  the  .Vth-order  predictor  except 
for  the  last  section.  This  feature  is  very  useful  in  handling  the  problem  of 
model-order  determination  and  reduced-order  modelling) ;  (iii)  a  cascaded 
structure  of  the  lattice  filter  (consisting  of  identical  sections)  which  is  very  con¬ 
venient  for  implementation  using  special  purpose  hardware,  microprocessors 
or  LSI  :  (iv)  in  normalized  versions  of  the  lattice  filter  ail  the  variables  are 
automatically  scaled,  making  it  possible  to  use  fixed-point  computations. 
(However,  normalization  sometimes  has  an  adverse  effect  on  the  numerical 
behaviour  of  the  algorithm  (see  Samson  and  Reddy  1982).) 

While  square-root  techniques  are  sometimes  applied  to  system  identification 
(Strejc  1980),  lattice  structures  are  apparently  not  used.  One  possible  reason 
is  that  efficient  recursive  algorithms  for  estimating  lattice  parameters  were 
developed  only  recently.  Another  reason  is  that  earlier  work  on  lattice  forms 
was  limited  to  all-pole  models,  while  most  realistic  plants  have  both  poles  and 
zeros.  The  work  (Morf  et  al.  1977,  I^ee  1980,  Lee  el  al.  1981,  Fried  lander  1982. 
Porat  et  al.  1981)  on  recursive  lattice  forms  provides  an  elegant  solution  to  the 
lattice  modelling  problem  for  both  all -pole  and  pole-zero  plants.  This  develop¬ 
ment  should  encourage  the  use  of  lattice  forms  in  system  identification  and 
adaptive  control. 

The  purpose  of  this  paper  is  to  present  lattice  implementations  of  the  follow¬ 
ing  commonly  used  recursive  parameter-estimation  algorithms :  recursive 
least-squares  (RLS).  recursive  instrumental  variables  (RIV).  extended  least- 
squares  (ELS),  and  recursive  maximum-likelihood  (RML)  (Xoderstrftm  et  al. 


g 
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1978,  Goodwin  and  Payne  1977).  The  idea  of  using  recursive  lattice  forms  for 
system  identification  was  previously  proposed  in  Morf  el  al.  (1977).  However, 
the  lattice  implementations  of  these  four  algorithms  were  not  discussed  in  full 
detail.  In  particular,  the  normalized  lattice  RIV  and  the  lattice  recursions 
for  arbitrary  model  orders  are  believed  to  be  presented  here  for  the  first  time. 

The  structure  of  the  paper  is  as  follows.  In  each  section  we  present  one 
of  the  lattice  algorithms  and  discuss  its  properties.  Owing  to  space  limitations, 
only  brief  derivations  are  included.  These  derivations  assume  some  familiarity 
with  the  projection  framework  for  developing  recursive  lattice  forms  which  is 
described  in  greater  detail  in  Lee  (1980),  Lee  el  al.  (1981),  Friedlander  (1982) 
and  Porat  el  al.  (1981).  We  have  attempted  to  make  this  paper  self-contained, 
but  some  of  the  background  material  has  been  deferred  to  the  references.  The 
results  in  this  paper  are  presented  in  the  normalized  case  only.  Unnormalized 
versions  of  these  algorithms  can  be  similarly  derived. 


2.  The  lattice  recursive  least-squares  algorithm 

In  this  section  we  consider  models  of  the  type  depicted  in  (I).  Given  a 
set  of  measurements  {»/0,  yT)  we  can  write 

X.VA.XM.  T^T  —  Uo:T  P  ) 


where 


n  0 


V  A 


L.y  r-l  •••  //  T-m  M  r- 1  •••  "  r-mj 
0r»the  estimate  of  the  parameter  vector  6  (eq.i.  (2)) 

.Vn:T~li/n<  I/tV 

The  least-squares  estimate  of  the  parameter  vector  is  given  bv  (indices  are 
omitted  for  notational  convenience) 

0TA(A"A')-'A"//,:T  (8) 

The  associated  error  vector  is  given  by 

*e:T  "  X6t  - 1/  -  A'(A"A*)-'A"  |y#;T  (9) 


where 


*o:T'“l*n . «r  I 


The  last  entry  «r  of  the  error  vector  is  given  by 

«T*w'f/  -  A'(A"A')-,A'  |//t.r 


where 


n  »|0 . 0,  II' 


Xote  that  Px  •  A'(A"JV)**A"  is  a  projection  operator  on  the  space  spanned  by 
the  columns  of  A',  i.e.  P  XP  x  —  Px-  In  Lee  (1980).  Lee  ft  al.  (1981),  Fried¬ 
lander  (1982),  and  Porat  H  al.  (1981 )  it  was  shown  that  projection  operators  of 
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this  type  can  be  recursively  updated  as  the  projection  space  X  is  changed  by 
the  addition  of  columns  x.  More  specifically,  the  following  update  formula 
can  be  derived 


U‘P*x+*  y  -  U'I*X  V  -  WPX*  x[x'/V  x]-‘x'PA-«  V  (11) 

where  V,  x  are  vectors  (or  matrices)  of  compatible  dimensions  and  P*  A  /  -  P. 
By  proper  choices  of  the  projection  space  X  and  the  vectors  V,  x,  the  un¬ 
normalized  lattice  recursions  are  obtained  (Lee  et  al.  1981,  Friedlander  1982). 
A  normalized  lattice  form  is  similarly  derived  by  considering  normalized 
projections 

p.y(F,  V)±[cpx*  r\-i‘*r'px*  v\  vpx*  V]~TI*  d2> 


which  obey  the  following  update  formula 


Px+AU,  F)  -  Px«\  x)p.v(*>  U)]-U'[px(<f'  V) 

-P.x(l\x)f>x(*'  F) ][/ — px(  V ,  x)p  A-(x,  F)|~™ 


(13) 


To  avoid  repeating  this  expression  we  will  find  it  convenient  to  define  the 
functions 


F(u ,  v,  tr) &[]-  irir’ ]-»'*[«  -  «•»][/  -  v'v)-Tli 
A’-1(«,  »,  ir)  A[/  —  mv'Y^uy  —  v'v]Tlt+tcv 


(U) 


Using  the  update  formula  (13)  we  can  derive  a  large  number  of  lattice  algorithms 
by  proper  choices  of  X,  x,  U,  V  and  proper  definitions  of  variables.  Table  I 
summarizes  the  variables  involved  in  the  LATTICE  RLS  algorithm.  The 
quantities  x„.r  and  x0.r"**  appearing  in  the  table  are  defined  by 


~y 

o  M 

0 

o 

o 

1 _ 

'  0  ’ 

0  n 

0 

X0:T  =* 

r  w.n  _ 

.  x0:T  * 

•  y»:Tm  * 

•j’o  «'o 

'Jo 

y 

j*  u 

r. 

Jf  T-m  T-m  . 

1 

* 

1 _ 

Ps(U,  V) 

S 

V 

V 

Comments 

*%.T 

Xp,o.T 

y«:T 

n 

«Yr 

^p.O.T-t-1 

uo:t 

n 

Joint-process 

lattice 

r*p.  T- 1 

Xp,O.T 

,p _ M 

Po:t 

n 

*V«.r 

Ap.o .T 

'Jo-T 

y,:r'+I 

*V,.r 

^p.»r+» 

U0’-T 

•Jo:t’*1 

*P.T 

Xm.p.T 

*o:r 

it 

Two-channel 

r».r-« 

*m.p.T 

*o:r 

n 

lattice 

^p+l.  T 

Xm.p.T 

X0‘-T 

n 

m&NA-XB+p 

Tabic  1.  Definitions  of  lattice  RLS  variables. 
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As  will  be  shown  shortly,  these  variables  define  the  lattice  structure  depicted 
in  Fig.  1.  where  0(A”)  is  an  operator  acting  on  the  input  vector  [«'.  r'|.  defined  by 
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Figure  I .  The  RI«S  lattice  form  :  (a)  overall  lattice  structure.  (A)  joint-process  lattice. 

(c)  two-channel  lattice. 


These  lattice  recursions  are  obtained  by  making  the  substitutions  in  (13) 
as  depicted  in  Table  2. 

The  following  facts  are  useful  in  interpreting  the  entries  of  this  table  : 


(i)  order  update 
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(ii)  order  and  time  update 
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Table  2.  Derivation  of  the  lattice  RLS. 


Another  point  that  requires  some  attention  is  the  interface  between  the 
joint-process  lattice  and  the  two-channei  lattice  in  Fig.  I.  Note  that 

f0.T^  PXpa- np.a.r^O-.T' ) 

r  (20) 

r0.T— 1  —  P.Y**-»«.o.r(X0.r''J-Aa+1’1-  flr)  ] 

If  we  assume  that  the  projection  operator  projects  first  on  ys  and  then  on  as 
(i.e.  the  square  roots  in  the  update  formula  and  the  definition  of  p  are  all  lower 
triangular),  then  we  note  the  following  : 

(i)  projecting  x0:T  on  Xxa_xiiot  involves  the  projection  of  n„.r  on 

(ii)  projecting  x0:Tx'i~v,,+l't  on  A'v<_vwor  involves  the  projection  of 
“o:r  on  A'.v^-.vw+i.r- 

Therefore  we  can  conclude  that 

«<l.  T  ”  I*I,/.V.4-.V«.  r*  *r'xA-XH.  T  r  I 

J  (21) 

rO.T-l  =  XA-XH+l.T-l  I  J 

We  can  now  summarize  the  lattice  RLN  algorithm  by  reading  off  the  proper 
entries  of  Table  2. 


Lattice.  RLS  algorithm 

The  algorithm  is  presented  here  for  the  case  XA&XB.  If  XA<XH. 
simply  interchange  y  and  u  and  XA  and  XB. 
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Input  parameters 

X  A ,  XB  —  model  orders 

A  -  exponential  weighting  factor 

<7  -  prior  covariance 
yT,  «  r  —  data  sequence 

Variables 

RT",  Rt*  *  estimated  eonvariance  of  y ,  « 
^%+i.T’  l.  r  “  reflection  coefficients 

reflection  coefficients 

«»p  T,  «p  r  »  forward  prediction  errors 

(dim  {y},  dim  {u},  dim  {[y',  w'J'}) 

'■"/#.  r-  r„.  t  ”  backward  prediction  errors 
(dim  {y},  dim  {(y',  «'  ]’}) 
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Initialization 

R_x***ol,  f?_,a=*o/ 

A\.-, P-1 . AM-AB 

A’*Pi.|*0,  p=l,  SA  —  XB+  1 
A'P.-i*r»..-i,“0'P=*1'  •••*  *Vfi 

.l/«in  loop 

At  each  time  step  do  the  following  : 

(i)  set 

Rt*  =*  AAr.i'  +  yrJf  7 
Rf* 15  -f-  MjtM  ^ 

«ar'  =  ^.T,'-(^Tv)-,/,yT 
«o.r“*(  Rt*)~  l/tur 
(i»>  update  joint-process  lattice 
For  p  =*  0,  min  {.Vri  —XB,  T) 

hl'l,+i,T!* r "  p.T- 1>  *Vr) 

rtn.T-i< 

r  “  ^Ir* t>' 7>_|,  **!,' t<  ^ 

r*  p<r,  «%,r) 

r>  r*i>.T'  A",,+i.rl 


omit  this  for 
p  =  min  {A’  A  -  A’ /?,  T*} 
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(iii)  set 

«#.r“  I**  XA-XH.T’  *"  t  I 

r0.r  =  (r*  XA-XIt. T‘  «“  A-J-A7i+».rl 

(iv)  update  two-channel  lattice 

For  p  =  o  to  min  {A'fl,  T  -  XA  +  XB)  -  I 

^V-l .T—  r  h.T- |t 

eM*i.r=  B(e».T’  ri>.T-i> 
r„~i.r**  F(ri,.T-i<  K'p+ur) 


Remark# 

As  the  lattice  recursions  are  started  they  may  involve  division  by  zero. 
It  can  be  shown  that  the  proper  procedure  is  to  set  to  zero  the  result  of  such 
division  in  the  scalar  case,  or  to  use  pseudo-inverses  in  the  matrix  case.  (See 
Porat  ft  al.  (1981)  for  details.)  For  coding  purposes  it  is  convenient  to  make 
/’’(•,  • ,  • )  and  F~l(  • ,  • ,  • )  into  subroutine  calls.  The  lattice  algorithm  then 
consists  of  repeated  calls  of  these  subroutines. 

Running  the  lattice  RLS  on  data  will  provide  a  set  of  reflection  coefficients 
parametrizing  the  plant  transfer  function.  In  some  applications  it  may  be 
desired  to  recover  the  estimates  of  the  {ait  b{)  parameters,  rather  than  to  con¬ 
tinue  with  a  lattice  structure.  Assuming  that  the  lattice  parameters  have 
converged,  this  can  be  done  by  looking  at  the  impulse  response  of  the  filter 
depicted  in  Fig.  1 .  Recall  that  this  filter  computes  the  normalized  prediction 
error  sequence  t,  where 

,-0  .-0 


where  [a(,  6(}  are  the  normalized  versions  of  {a(,  b(}.  Xote  that  the  impulse 

50,  ff, . SXA)  while  the 

The  unnormalized  para- 


response  from  the  y  input  to  the  < 

outgut 

will  be 

impulse  response  from  «  to  «  will  be 

(V  K  ■ 

...  5V«} 

meters  can  be  obtained  by  setting 

«(  =  «(Tl«o 

» *  * .  ■ 

...  XA 

b{  =  aB~l  b(. 

*  *  o. . 

...  SB 

(23) 


This  method  of  computing  {a*,  b( }  from  the  reflection  coefficients  does  not 
give  the  exact  least-squares  estimates  of  these  parameters.  However,  if  the 
reflection  coefficients  have  converged  ‘  sufficiently  ’  these  estimates  will  be 
very  close  to  the  optimal  estimates.  A  slightly  more  complicated  lattice 
filter  is  available  for  computing  the  exact  least -squares  estimates  of  {a,-,  b{), 
from  the  information  provided  by  the  lattice  RLS  ;  see  Friedlander  (1982) 
and  Porat  el  al.  (19811  for  a  more  detailed  discussion. 

The  lattice  RLS  differs  from  the  standard  RLS  algorithm  in  several 
respects  . 


(i)  Initialization.  The  order  recursive  nature  of- the  lattice  RLS  makes  it 
possible  to  eliminate  transient  phenomena  caused  by  incorrect  initial  conditions, 
leading  to  faster  startup. 


•ovievvy.’ 


-  -  ’  V  ' 

1  - 
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(ii)  Normalization.  All  quantities  the  lattice  RLS  (with  the  exception  of 
R*r,  Rmt)  have  magnitudes  less  than  one. 

(iii)  Computational  requirements.  The  lattice  RLS  is  computationally 
efficient,  requiring  O(X)  operations  per  time  step,  where  X  =  XA  +  XB.  The 
usual  implementation  of  the  RLS  requires  0(A’*)  operations  per  time  step. 
Since  square-roots  are  time-consuming  operations  on  general-purpose  com¬ 
puters,  the  efficiency  of  the  lattice  forms  becomes  apparent  only  for  fairly  large 
values  of  X.  However,  implementations  on  special-purpose  hardware  designed 
to  take  advantage  of  the  lattice  structure,  can  be  very  efficient.  Note  that 
the  computation  of  the  {at,  b()  parameters  from  the  reflection  coefficients 
requires  0(A'*)  operations.  This  computation  can  be  avoided,  however,  by 
reformulating  the  problem  for  which  the  parameters  were  estimated  so  that  it 
it  will  use  directly  the  reflection  coefficients. 

(iv)  Order-recursive.  The  lattice  RLS  is  not  only  time-recursive,  but  is 
also  order-recursive.  This  makes  it  possible  initially  to  overdetermine  the 
plant  order  and  to  choose  a  lower-order  model  after  the  parameter  estimates  are 
computed. 

Finally  we  note  that  the  lattice  algorithm  presented  in  this  section  was 
only  one  of  many  different  lattice  forms  (the  so-called  normalized  pre-windowed 
form).  The  unnormalized  lattice  recursion  and  the  convariance  lattice  form 
are  presented  in  Lee  el  al.  (1981)  and  Porat  et  al.  (1981).  The  pre-windowed 
form  is  simpler  than  the  covariance  form  and  is  probably  better  suited  for 
system  identification  applications. 


3.  The  lattice  recursive  instrumental  variable  algorithm 

The  RLS  algorithm  provides  biased  estimates  when  the  disturbance  process 
e,  (see  eqn.  (1))  is  non-white.  The  intrumental  variables  method  (Young  1970. 
SoderstrOm  and  Stoica  1981,  Wong  and  Polak  1967)  was  derived  to  eliminate 
this  problem.  The  parameter  estimate  dT  is  given  by 

$T-(Z'X)-'Z'g0:T  (24) 

where  Z  is  an  instrumental  variable  matrix 


Z 


m,  it,  T ' 


0  0 


y  o 


y » 


y  T- 1  y  T-m  u  T-l  “  r-«i  J 


(25) 


The  instrumental  variables  yt,  u,  can  be  chosen  in  different  ways.  A 
typical  choice  is  to  set  £,=  u,  and  y,  to  be  the  output  of  a  filter  driven  by  u„ 
for  example 

xa  xu 

//«=-  I  «<.?/-,+  I 

•  •  l  *  -  0 


(26) 
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Note  that  the  projection  formula  (10)  is  replaced  in  this  case  by 

(27) 

which  involves  the  non -symmetric  ‘  projection  operator '  pzx.  By  analog}* 
with  the  derivation  outlined  in  §  3,  we  define  a  normalized  projection  operator 

PzXl(’.V)HL:'PzxlU)-ll9l'Pzx°V\VPzx'V]-™  (28) 

In  Appendix  A  we  prove  the  following  update  formula  for  this  operator 
»  V)  =  \!  ~  Pzx(t  ' <  2)Pz\ ~ z)Pzx(X'  ^ ‘ 

x  IPz.vl^  >  ~  Pzx(^  <  z)Pzx~l(x<  z)Pzx(X’ 

X  [1  -  Pzx(  V-  2)pZX~l(X,  Z)pzx(X,  V)]-™ 

To  avoid  repeating  this  complicated  expression  we  will  find  it  convenient  to 
define 

F{u,  v.  ir,q,  r,  *)  =  |/  —  qa'yr\-xlt\u  —  —  w*«-1t>|_r/*  (30) 

and  its  ‘  inverse  ‘ 

F~l(v,  v,  ir,  q,  r ,  «)»[/ —  ga-,r|l/*M[/  —  «•«-*»  |r/*  +  g«"lr 

Using  the  update  formula  (29)  we  can  derive  several  versions  of  the  lattice 
RIV  by  proper  choices  of  Z,  X.  z.  x,  C.  V.  To  simplify  the  presentation  we 
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consider  in  this  paper  only  the  case  a  here  XA  —  SB.  The  more  general  case 
can  be  derived  in  a  straightforward  manner  following  the  steps  outlined  in 
§  2.  Table  3  summarizes  the  variables  involved  in  the  recursions.  We  will 
also  use  the  following  definitions 
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The  recursions  are  obtained  by  making  the  substitutions  depicted  in  Table  4, 
in  the  update  formula  (29). 

The  following  facts  are  useful  in  interpreting  the  entries  of  this  table  : 


(i)  order  updates 

X  jj+l.T  *  -V  r  +  T  *  %II.  T  +  *o :TI>+1 

(ii)  time  and  order  updates 

X,,* 1.  T+l  =“  X T  +  X'o-.T<  Z/j+l.  r-l  =  %I>.  T  +  Z*:T 

(iii)  lime  update 

f". 


o  .  o 


Note  also  that  pZA-(f\  F)=*p' vz(  F, 

We  can  now  summarize  the  lattice  RIV  algorithm  by  reading  off  the  proper 
entries  of  Table  4. 


Lattice  Rl  V  algorithm 

We  denote  here  M  =  dim  {(y'r,  «'r|' j. 


Input  parameters 

S  m  model  order 


A  *  exponential  weighting  factor 
a  m  prior  covariance 
i/T,  Nr*data  sequence 
f/r,  uTm  instrumental  variable  sequence 


■  A-‘ -s  •  ^  ■  •- 


>  >*•  ,•  *  vvS» 

„  -  V  \  or  0  o  O  m3. 
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Variables 

Rt ^r'-estimated  covariance  of  [y’r.  [y'r»  «'r]'  x  J/) 

A’xi,  Alj,  Kxt,  Klx,  Klt,  A'**,  Ax*,  A.'11  =  reflection  coefficients  (J/  x  J/) 
^p.r-  ^/..r*  r*//,r*  r  *  prediction  errors  (M  x  1) 

«%.r*  r*  ?Vr,“aux*l**ry  prediction  errors  (M  x  1) 

Initialization 

<r/ . 

Reflection  coefficients  and  backward  prediction  errors  are  all  initialized 
to  zero. 


loop 

At  each  time  step  do  the  following  : 

(i)  Set 

RT*mhRT-ix+\y'T,  «'ri ¥r.  «VI 
V-  «'rl'[fr,  S'rJ 

‘t/-'’zrz-^rx-\r,-lV)*,,Vr.  «V1' 

«ar’-''o.r,-«,.r1  - \ts •  (RT,)~l/t[y' T.  *'tY 
(ii)  For  p-0,  ....  min  {A*.  T}  —  1,  do  : 

A’*Vur“  '“/..r-i.  «Vr.  «%.r.  7) 

^Vi.r“^*,(^Vi.r-i'  '“p.r-i-  rVr-i>  «%.r>  0 


A'xVi.r“ 

^-'(A'xV,r- 

-I* 

pp.t- 

■V  r*p,T-l>  ** p,T'  *P. TX<  l) 

A’l*i/+i.r=s 

J,-1(A',*#<+I  r, 

C3 

/<.  r-  ^ 

p.t*  rtp,r-i>  ^.r-i'  ^ ) 

A'!* 

11  1/4- 1,  r 

^-‘(A'lV..r- 

*  p,r> 

^e.T’  r*i/,r-i>  *%.r-i<  *) 

I'i: 

A  p+i,r 

-I» 

^p.T’  ^n, T~i<  **p, r-i.  0 

A’r:  ™ 

n  /j+I*  r 

^,-,(A"p4-i  ,r- 

1* 

**p.r> 

*"*//,  T’  *Zp,T'  I) 

71  p^i.r2* 

*Vr- 

i>  r-i>  r-i>  r-i* 

/) 

** 11,  T  - 

-I* 

^,/.r 

-if  ^*Vn,  r»  ^*Vn,r» 

r%-i.r” 

r» 

?  p.T< 

1  ^  p+i.r*  ^ 

*'(/-!.  r* 

r ’  ^//,  t*- 

1* 

r„.T- 

-1*  K*s„+i,T'  a** 

//4-t,r) 

^(^/..r-n  rrp.r<  r  p.r>  i,r>  A*x%+i,r) 

*"*  /i.r-i>  A'ix  ,,+i.y  A**s  p+i.r) 

r*ii+l.T  m  *”/»,  T’  **  p,T'  Ri*  |/4.|,  r-  Ai!  j/4-l.r*  Kx:  n+i'T) 

*Vi.r’^(r»,r’  **  u.t~ i>  A"  A'5* , 
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As  can  be  seen  from  these  equations,  the  normalized  lattice  R1V  is  fairly 
complex.  This  is  due  to  the  more  complicated  projection  update  formula  and 
to  the  fact  that  various  identities  which  were  true  for  the  symmetric  projection 
operator,  no  longer  hold  (e.g.  p2A(t',  V)  +  p'MS(V,  l')).  The  unnormalized 
version  of  the  lattice  RIVr  turns  out  to  be  considerably  simpler  than  the 
normalized  recursions  presented  above  (see  Appendix  B).  This  is  different 
from  the  situation  in  the  lattice  RLN  where  the  normalized  version  is  the  simpler 
one.  The  {a(,  6{}  parameters  can  be  recovered  by  computing  the  impulse 
response  of  the  variance  normalized  lattice  RiV. 

Finally  we  note  that  the  unnormalized  lattice  RIV  has  been  developed 
independently  by  several  authors  (see,  for  example,  Samson  1982,  Cadzow  and 
Moses  1981).  An  approximate  lattice  RfV  was  presented  by  Prevoeto  et  al. 
(1982). 

4.  The  lattice  extended  least-squares  algorithm 

In  this  section  we  consider  the  following  ARMAX  model 

.vj  .vs  .vc 

Jfi«-  I  I  Mi-<+  I  c(r,_,  +  e,  (31) 

■  -  I  0  im  i 

where  o,  is  an  unmeasurable  white-noise  disturbance  process.  If  it  were 
possible  to  measure  v„  this  would  have  been  a  standard  linear  regression 
problem,  and  the  RLS  algorithm  could  be  applied.  The  ELS  method  is  based 
on  the  idea  of  replacing  v,  by  its  estimate,  the  prediction  error  c,  (SOderstrOm  et  al. 
1978,  Panuska  1989,  Solo  1979).  The  lattice  ELS  will,  therefore,  consist  of 
two  steps  :  (i)  use  the  lattice  form  as  a  prediction  filter  to  compute  *, ;  (ii)  use 
the  lattice  RLS  for  the  known  input  case  (with  y„  u„  »,-«,)  to  update  the 
parameter  estimates. 

To  describe  the  lattice  ELS  for  the  model  presented  above  we  must  first 
present  the  basic  update  formula  of  the  lattice  RLS  for  ARMAX  models. 

Lattice  RLS  (ARMAX)  algorithm 

We  assume  that  X A  2  XB  ^  XC.  For  other  cases  we  simply  have  to  reorder 
the  inputs  y,,  u„v,  so  that  the  corresponding  model  orders  appear  in  decreasing 
order. 

A'*.  A'“,  A’r,  R,  Rr,  K  «  reflection  coefficients 

«*.  «r.  «,  ir,  ««■  forward  prediction  errors  (dim  {y},  dim  {«},  dim  (»}. 

dim  {(y',  «’]},  dim  {r},  dim  {[y\  »'1» 

r»,  r.  r  *»  backward  prediction  errors  (dim  {y},  dim  {[y\  «' |). 
dim  {([y',  »']}) 

(i)  for  p»0 . X A  -  XB  (or  up  to  T,  during  start-up) 

R~l(R*r-ri.T-i<  r *  „,t- i>  «%.r) 

i,r“  fVr- 1* 


skip  for  p  «  XA  -  X  B 
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" *.r>  «Vr> 

^Vi.r,*^'l(^ViT-i> r,f  ».*••  «%.t> 

«V«.r- /'(•Vr.  rV r>  A'Vi.r) 

‘Vi. r  *  A(*%,  j*.  ^V, j>t  ^Vi r) 

(ii)  set 

.Y.O-.VJ*  *"  .v.4-.v«+i.r]" 

?#.r“[rl'  xa-sh. r»  .v^-.v/i+url* 

***.  p  "  ‘V-.v/zti.  r 

(iii)  for  p  •  0  to  .VB- A’C  (or  up  to  T  -  A’ ^  +  A’ if  during  start-up) 

Rl>+l.Tm  F  l(Rp+l.T-l>  ^  p.T- 1*  *|#,f) 

‘Vtr-^lVr-i.  R„+ut)  skip  for  p- NB~  AT 

V*i.r"  ^V*ur) 

^V+i.r*  ^-,(^Vn.r-i*  ^  «.r>  ** p.r) 

*Vn.r®  A(#* PiT,  pSp.y,  ^%+i.r) 

(iv)  set 

*«.  r  **  l*  «,  r-  **  .vw-.vc*i  I* 

Vr" l^’o. r*  **  NB—xc+i  1* 

(v)  for  p  -  0  to  A’C- 1  (or  up  to  T-AM+A’C-  I  during  start-up) 

M^+i.r-11  r p,r-t»  Vsr) 
ti>*i.Tm  F(*i,,T’  V.r-i*  A|,+i.r) 
r#.+i.r®  «p.r>  K'M+l  T) 

The  corresponding  lattice  structure  is  depicted  in  Fig.  2.  The  detailed  structure 
of  the  various  sections  is  similar  to  that  of  the  lattice  RLS  of  Fig.  1,  with  some 
obvious  modifications. 


(SCI 


Joint  Proem  Two-Ctunml  Tlir*e-Ch*nn*l 

Lottie*  Joint  Proem  Lottie* 

Uttlc* 


Figure  2.  The  ARMAX  RLS  lattice. 
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The  lattice  ELS  algorithm  can  now  be  summarized.  The  RLS  lattice 
update  described  above  will  be  used  in  two  modes.  One  is  a  prediction  mode, 
in  which  all  the  reflection  coefficient  updates  are  skipped  and  the  old  valuesof 
the  state  variables  {r»p.r_ t,  r* ^.r-i*  t»,t- il  a*®  retained.  The 

second  is  a  regular  update  mode  in  which  both  reflection  coefficients  and  state 
variables  are  updated. 

Lattice  ELN  algorithm 

All  matrix  square-roots  are  lower  triangular. 

I  initialization 

All  reflection  coefficients  and  backward  prediction  errors  are  initialized 
to  zero. 


Main  loop 

At  each  time  step  do  the  following  : 

(i)  compute  prediction  errors 
A/fT_i*'  +  yry'r 
/?j.“  *  +  *r**  t 

<0,T*~(*TU)-Vt*T 


Call  lattice  RLS  update  (ARM AX),  in  prediction  mode 


(ii)  set 

tir-  last  entry  of  «  Vr.r 

(iii)  update  lattice  variables 


+  VfV  f 

')'vt6r 

Call  lattice  RLS  (ARMAX)  in  update  mode. 

As  before,  the  parameter  estimates  {a(,  h(,  c(\  can  be  recovered  by  looking  at 
the  impulse  response  of  the  lattice  form  depicted  in  big.  2,  from  the  inputs 
y,  n,  v  to  the  prediction  error  output.  For  a  more  detailed  discussion  of  the 
lattice  ELS  for  the  case  of  ARMA  (.V,  .V)  models  (i.e.  A  A* AC.  A  R»0),  see 
I<ee  el  al.  (1980).  An  approximate  lattice  ELS  algorithm  for  general  ARMA 
processes  was  presented  in  Benveniste  and  Chaure  (1981). 


5.  The  lattice  recursive  maximum-likelihood  algorithm 

The  RML  algorithm  has  improved  asymptotic  convergence  properties 
compared  with  the  ELS  algorithm  presented  above.  To  simplify  the 
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presentation  we  consider  in  this  section  ARMA  models  rather  than  the  more 
general  ARMAX  models,  i.e. 

.V.l  .VC 

»#--  Z  aiUt-i+  Z  '•(«'»-< +  (32) 

.-i  .-i 

where  v,  is  a  white-noise  process.  The  RML  algorithm  (SOderstrOm  et  at.  1978. 
SoderstrOm  1973.  Ljung  1977,  1981)  can  be  summarized  as  follows.  Let 

9  » (a, . a  v  ,.  r, . rxr  jT  *  parameter  vector 

4i-l  . i.  '#-i . e,_.vrlT“d*ta  vector 

.  -f/i-x.t-  *i-i . filtered  data  vector 

where  e,.  e,.  ij,  will  be  defined  later.  The  update  equations  are 

6,-6t-K  + 

(33) 

*V* -/•»„■ •  +  r 

The  filtered  quantities  are  obtained  by 


y, «l !/<?,(*)  ly i  J 

where 

<?<(z)«l+*,(<)2-'+  ...  +<V(.(t)r-vr  (35) 

where  z~y  is  the  unit  delay  operator.  Note  that  (33)  can  be  rewritten  as 

fV'  6,-Pc1  (36) 

Let  us  define 

^«-i  (37) 

and  sum  up  the  difference  equation  (39)  to  get 

PCl6,~  Z  0c r,  (38) 

f-l 


^■f  Z  M<T  1  Z 

L  <-i  J 


Equation  (39)  can  be  recognized  as  the  solution  to  the  problem  of  estimating 
the  process  .r,  from  the  components  of  the  vector  </>,.  Thus,  given  the  variables 
we  can  apply  the  recursive  least-squares  algorithm  to  estimate  the 
parameter  vector  9.  The  joint  estimation  problem  described  above  can  be 
solved  by  the  following  joint-process  lattice  form. 

Joint-process  tiro-channel  lattice 
We  assume  here  XA  -  AY'-  X. 
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Main  loop 

R?  m  •+■  ZfZ  f  R^  m  +  A  (40  a) 

*%Tmiro.Tm  Rr~ut  %t  Ur'-lVr1''^  (406) 

For  j»*l)  to  «\\  do 

^r*t.T~  F  *(^j*4’1.t— i>  ^  m. r— I*  *p, r)  (41**) 

^Vi.r*^'l(^Vi.r-i>  r  *>.r>  «* *.r)  (41 

V*i.r*  rD.r~ i*  K,+ut)  (41  c) 

F(ri*.T-i’  V r*  (4*  <0 


The  algorithm  described  above  will  both  update  the  parameter  estimates 
(reflection  coefficients  A',  A-*)  and  filter  incoming  data  to  compute  a  set  of 
prediction  errors.  Sometimes  we  want  to  use  this  lattice  structure  for  filtering 
only.  In  this  case  only  (40  6),  (41  c),  (41  d )  and  (41  e)  need  to  be  used.  To 
distinguish  between  these  two  cases  we  will  call  LA  TUP  the  algorithm  that 
performs  the  full  computation  (40)-(41)  and  LATFIL  the  algorithm  that  does 
filtering  only. 

The  lattice  structure  described  above  can  now  be  used  to  implement  the 
RML  algorithm  (Friedlander  et  at.  1981).  This  will  require  several  steps  of 
filtering  and  parameter  updating.  The  following  set  of  equations  summarizes 


the  RML  algorithm  in  non-lattice  form 

«r-.Vr+  I  *<(T-nyr-,-  (42*) 

<■1  (•! 

-fT.0+  £  djT- l)jT-i-  Z  •JilT’-IXr-,  (426) 

*  •  I  #  •  I 

xT»«r  +  iT  ;  perform  least-squares  parameter  update  (42  r) 
*TmStT+'  Z  J)  (42 rf) 

»  •  I  «  •  I 

.VC 

yT-yr  +  o-  £  tATtfr.i  (42  e.) 

i  ■  I 

<T-*T  +  »-  £  f({T)iT  i  (42/) 

»■! 


fc 


i 

•v"’ 

V' 
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Call 

Y 

V 

X 

Equation 

LATFIL  (1) 

y»-\ 

*1- 1 

yt 

(42  a) 

LATFIL  (2) 

#1-1 

3<-l 

0 

(42  6) 

LATUP  (2) 

9t- 1 

*1-1 

(42  c) 

LATFIL  (1) 

Yi-i 

*1-1 

y, 

(42  d) 

LATFIL  (3) 

0 

#1-1 

y» 

yt 

(42  e) 

LATFIL  (4) 

0 

1 

« / 

it 

(42  /) 

Table  5.  The  RML  lattice. 


This  set  of  equations  can  be  implemented  by  repeated  calls  of  the  lattice 
form  described  above.  The  input  and  output  for  each  lattice  call  are  summarized 
in  Table  5.  Note  that  four  different  ‘  state  vectors  ’  need  to  be  stored  corre¬ 
sponding  to  0,  <l>  and  the  pre-filters  for  y  and  e.  These  four  cases  are 
distinguished  in  Table  a  by  the  index  of  the  lattice  call  (for  example,  LATFIL  (1 ) 
represents  the  filters  with  yT,  e.T  as  inputs,  while  LATFIL  (2)  has  yT,  iT  as 
inputs). 

6.  Conclusions 

The  lattice  equivalents  of  several  system  identification  algorithms  were 
presented.  These  algorithms  provide  a  computationally  efficient  recursive 
solution  of  linear  least-squares  estimation  problems.  In  the  area  of  digital 
signal  processing  lattice  filters  are  often  preferred  over  their  tapped-delay-line 
equivalents  because  of  their  relative  insensitivity  to  roundoff  errors.  In 
adaptive  processing  applications,  lattice  filters  have  shown  improved  con¬ 
vergence  behaviour  compared  to  the  popular  Widrow-Hoff  LAIS  algorithm, 
lattice  structures  also  lead  to  processing  architectures  that  are  quite  different 
from  those  related  to  the  RLS  and  similar  algorithms.  This  modular  pipe¬ 
lined  architecture  has  potential  advantages  in  hardware  and  VLSI  imple¬ 
mentations  of  the  algorithms. 

Relatively  little  work  has  been  done  in  the  application  of  lattice  forms  to 
system  identification  and  adaptive  control.  Considerably  more  analysis  and 
simulation  studies  are  needed  to  assess  the  usefulness  of  the  techniques  pre¬ 
sented  in  this  paper.  Of  special  interest  would  be  tests  performed  on  finite 
word  length  machines  and  plants  with  high  order  dynamics.  These  conditions 
often  lead  to  numerical  problems  in  standard  recursive  parameter  estimation 
algorithms.  It  is  hoped  that  this  paper  will  stimulate  research  in  this  area. 

Appendix  A. 

Derivation  of  the.  update  formula  for  non-sym  metric  projection  operators 
Definitions 

Ft,*  AZ(A"Z)-'A" 

P'z.xH-Pzx 
Z  +  zh\Z  z]  A'+rA|A'  r| 


V 

h  _ 
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Update  formula  for  PT 

P'z+t.X-x  *  zx  ~  P* xx  tfa'P* zx  *1  ~lx'Pr zx 

Proof 


x'z  jrvHrjr 

3Z*t,X+Tm(Z  2  I 

x'Z  x'z  J  x' 


Use  the  following  matrix  identity 

[A  B]->  [d-*  01  r-^-‘fTJ 

+  A0-M-<u-'/| 

CD  J  L  (»  L  /  J 

A  n±D-CA-'H 
to  invert  the  piatrix  in  (A  2)  to  get 

pz++x+*-WZ\~x  X'  +  (/  -  Z(X'Z)~xX')zX~lx'(I  -  Z(X’Z)-*X’)  | 

A  -  x'[l  —  Z(X'Z)~lX‘  ]2  j 

Equation  (A  1)  follows  directly  from  (A  4)  and  the  definitions. 


(A  I) 


(A  2) 


(A  3) 


(A  4) 


Xormalization 

p*a(<\  V)±[cp«zx  u\-iitu’p*zx  v\  rp*MX  v\-*<* 

\<"P*zx  t'l -tnr-P'z^x+*  v\V'P*zx  V]-Tlt 


Set  VmC 


” Pzx(t  ’  V)  —  Pzx(l! ■  z)p~xzx(x<  2)Pzx(x<  V) 


[rp^zx  f'i-,/,r  P‘z*t..Y*x  r\rp*tx  I-}-™ 

mI  —  Pzxil  •  z)Pzx~l(2'  x)Pz. v(x»  ^  ) 

Set  V  -  r 

W'P*ZS  VY'nrP*s+*x+*  V\VP*zx  v\-Tlt 

=  /  —  Pzxl  ^  z)pzx~l(z’  *)Pz.v(x> 

Normalizing  CP^zx  ^  hy  the  square-roots  of  the  last  two  equations  gives 
pz^x+av,  v)-\i- pzx(i\ z)pzx~'&> z)pzx(x>  w* 

x  Ipzxi l  >  V)-PzxW.  z)PgX-H*>  z)Pzx(x>  V)\ 

x\!-  Pzx( V'  *)Pzx~l(*'  z)Pzx(x-  V)~Tlt 


Appendix  B 

The  unnormalized  lattice  R1  V 

The  unnormalized  lattice  R1V  consists  of  two  RLS-type  lattice  filters,  as 
depicted  in  Fig.  3.  The  input  to  the  upper  lattice  is  the  data  sequence  x ,  and 
to  the  lower  lattice  the  instrumental  variable  sequence  z,.  The  reflection 
coefficients  of  the  two  filters  are  determined  by  a  common  set  of  coefficients  : 
A’*\  K*1,  Ku,  Kn.  Table  «  summarizes  the  definitions  of  all  the  variables  in 
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the  lattice  RIV.  Making  the  proper  substitutions  in  the  update  formula  for 
r'/*zv  V  gives  the  set  of  lattice  recursions.  The  necessary  substitutions  are 
summarized  in  Table  7.  Reading  off  the  entries  of  this  table  gives  the  following 
equations 


*X/»+i.ras*'c h.t~  ^ lr'/.,r-i 

.T“<c n.T~  ,,+  t i^l.T^'^H.T-l 


lattice  I 


lattice  2 


n+i  .r  =  ^/,.r 


r=  +rr/».r-i  r"  /(.r-i/l*  — V/i-i.r-i) 

^x:i,-i.tss  Mtxx,,+i.T-i  +(Xi>.t  *'  —  y/<-i.r-i) 

y„.r=!y/.-i1r-i  +  «"  ) 


(B  I) 
(B  2) 
(B  3) 
(B  4) 
(B  5) 


(B  6) 


time  updates 


(B  7) 
(B  8) 


|  time  and  (B  9) 
>  order 

A’«;i+tr*t  =  A’Vl.r+^V1.r(A'Vl.r)-,A’‘-V.l.r  ]  update  (BIO) 


K xl„+t.  r  =  A  XI„+i.  r  +  ^  xij#+i.  T^x'e*i.  r)-1  ^x:e~i.  t 
y„.T-l  xy,,-X.T  +  r*  i,.T-l(Kx-„<.x,r)-'r'i..T~l 


order  update 


(B  II) 
(B  12) 


Xote  that  we  have  introduced  the  exponential  weighting  factor  A  into  the 
time  update  equations  (B  5)-(B  8). 

The  complete  RIV  algorithm  can  be  implemented  in  several  ways  using 
these  equations.  For  starting  up  the  algorithm  it  is  necessary  to  use  the  order 
(or  time  and  order)  update  equations  (B  10)  and  (B  12)  for  Kx:.  Ki:.  After¬ 
wards  the  time  update  equations  (BO)  and  (B  8)  will  be  used  instead.  The 
initial  conditions  for  the  algorithm  at  time  step  T  are 


o.  T ' 


n,  T  ’ 


o,  T 


0.  T  ■ 


y-i.r  = 


^  *'o.t  “  R -  £  rTz  t  —  ^  r— i  +  xtz  i 


Before  start-up  all  the  reflection  coefficients  and  state  variables  are  set  to 
zero.  A  complete  description  of  the  unnormalized  lattice  RLN  can  be  found  in 
Lee  (1980),  Lee  ft  al.  (1981)  and  Friedlander  (1982).  A  comparison  with 
cqns.  (B  I  )-(B  12)  leads  to  one  possible  implementation  of  the  lattice  RIV. 
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